import numpy as np
import matplotlib.pyplot as plt

Lx=0.1 # taille du domaine de travail en mètre
Ly=0.1

Nx=100  # nombre de points de la grille
Ny=100

dx=Lx/Nx # intervalle entre deux points
dy=Ly/Ny

precision = 0.5 # règle la précision du calcul : si la methode converge à ce
                 # nombre près, le programme arrête les itérations

V=np.zeros((Nx,Ny))       # tableaux à 2 dimensions vide au départ
V2=np.zeros((Nx,Ny))      # V : potentiel
B=np.zeros((Nx,Ny))       # Ex et Ey composantes du champ électrique
Ex=np.zeros((Nx,Ny))      # B : si B=1, le potentiel est imposé, B=0 vide
Ey=np.zeros((Nx,Ny))

# conditions aux limites : électrode à 5 kV dans une boîte passant en
# bas de 5 kV à 0V en haut

B[0,:]=1
B[:,0]=1
B[Ny-1,:]=1
B[:,Nx-1]=1

V[0,:]=5000

V[:,0]=np.linspace(5000,0,Ny)
V[:,Nx-1]=np.linspace(5000,0,Ny)

V[Ny-1,:]=0


Nx2=abs(Nx/2)

for i in range(Ny-1):         # électrode  = demi élipse portée à 5000 V
    for j in range(Nx-1):
        if i**2+(j-Nx2)**2*100<50**2:
            V[i,j]=5e3
            B[i,j]=1



#------------calcul du potentiel méthode de Gauss Seidel -------------------
emax=1000                       # on met une valeur arbitraire > précision
k=0                             # pour que le while ne stoppe pas tout de suite
                                # les itérations, emax sera l'écart max de
                                # potentiel entre deux itérations

w=2/(1+np.pi/max(Nx,Ny))        # paramètre optimal pour la convergence

while emax>precision:           # tant que la présicion n'est pas atteinte
                                # on met à zéro au départ pour calculer le max
    emax=0                      # d'écart sur l'ensemble de la grille

    for i in range(Ny):
        for j in range(Nx):


            if B[i,j]==0: # si on est dans le vide : on résout l'equ de Laplace

                V2=V[i,j]
                V[i,j]=(1.-w)*V[i,j]+w*(V[i+1,j]+V[i-1,j]+V[i,j+1]+V[i,j-1])/4.

                # on travaille directement sur V et on utilise les valeurs
                # de la grille en cours de calcul pour faire converger la
                # méthode plus vite
                # plus besoin de travailler sur une matrice copiée !

                emax=max(emax,abs(V2-V[i,j])) # on compare pour trouver
                                                   # l'écart maximum
    k=k+1                 # on compte le nombre d'itérations
#---------------------------------------------------------------------------

print(str(k)+" itérations")

# -------------calcul du champ électrique Ex et Ey--------------------------
Emax=0

for i in range(1,Nx-1):
    for j in range(1,Ny-1):

        Ey[i,j]=-(V[i+1,j]-V[i-1,j])/(2*dy) # calcul du gradient selon Oy
        Ex[i,j]=-(V[i,j+1]-V[i,j-1])/(2*dx) # selon Ox
        Emax=max(Emax,(Ex[i,j]**2+Ey[i,j]**2)**0.5)
#----------------------------------------------------------------------------

print("Emax  = "+str(round(Emax))+ " V/m" )

# tracé de la carte de champ : potentiel et champ électrique

x=np.linspace(-Lx/2,Lx/2,Nx)
y=np.linspace(-Ly/2,Ly/2,Ny)
X,Y=np.meshgrid(x,y)


plt.figure(figsize=(10,8))
cf=plt.contourf(X,Y,V,100,cmap='jet')
graph=plt.contour(X,Y,V,12,colors='black')
plt.clabel(graph,inline=1,fontsize=10,fmt='%3.2f') #affiche les valeurs de V
plt.colorbar(cf)
plt.xlabel("en m")
plt.ylabel("en m")
plt.title("Carte de champ et potentiel")
plt.streamplot(X,Y,Ex,Ey,density=1,color='black')
plt.show()

normeE=(Ex**2+Ey**2)**0.5
plt.figure(figsize=(10,8))
plt.plot(y,normeE[:,50])
plt.grid()
plt.xlabel("en m")
plt.ylabel("en V/m")
plt.title("Norme du champ électrique - coupe plan vertical")
plt.show()