### Resolution de l'equation de Laplace sur un carre tel que le condensateur est centre

## Importation des modules
import numpy as np
import matplotlib.pyplot as plt

## Parametres du probleme en USI
L = 50     # 1/2 cote du carre etudie en m 
L_cond = 50 # Longueur des armatures du condensateur
e_cond = 15 # Distance entre les deux armatures

# Valeur du potentiel sur les bords en V
V0 = 0
# Valeurs du potentiel sur les armatures en V
Vplus = 5
Vmoins = -5

## Discretisation de l'espace
# Nombre de cellules selon chaque direction et pas d'espace en USI
n = 100
h = L/n   
  
## Tableau des potentiels et CL
# Création du tableau des potentiels et du tableau de booléens représentant les frontières
V = np.zeros((n+1,n+1))
B = np.zeros((n+1,n+1), dtype=bool)
# Affectation des CL
# sur les bords
for i in range(n+1):
    V[0][i] = V0
    B[0][i] = False
    V[i][n] = V0
    B[i][n] = False
    V[n][i] = V0
    B[n][i] = False
    V[i][0] = V0
    B[i][0] = False
# sur les armatures
ihaut = (n-e_cond)//2
ibas = (n+e_cond)//2
jgauche= (n-L_cond)//2
jdroit = (n+L_cond)//2
V[ihaut][jgauche:jdroit] = Vmoins
B[ihaut][jgauche:jdroit] = True
V[ibas][jgauche:jdroit] = Vplus
B[ibas][jgauche:jdroit] = True

## Algorithme de Jacobi
def ecart(V1,V2):
    diff = 0
    for i in range(n):
        for j in range(n):
            diff += (V1[i][j] - V2[i][j])**2
    e = np.sqrt(diff/(n**2))
    return e

def iteration_jacobi():
    oldV = V.copy()
    for i in range(1,n):
        for j in range(1,n):
            if B[i][j] == 0:
                V[i][j] = (oldV[i+1][j] + oldV[i-1][j] + oldV[i][j+1] + oldV[i][j-1])/4
    e = ecart(V,oldV)
    return e

def jacobi(eps,imax):
    e = 2*eps
    i = 0   # compteur d'itérations
    while e > eps and i<imax:
        e = iteration_jacobi()
        i += 1
    return i

## Maillage de l'espace
# Listes des xi et des yj 
xlist = np.linspace(-L,L,n+1)
ylist = np.linspace(-L,L,n+1)
# Matrice correspondant au maillage de l'espace 2D
X, Y = np.meshgrid(xlist,ylist)

## Appel de la fonction pour resoudre et trace des equipotentielles
N = jacobi(1e-3,1e3)
print(N)
plt.contour(X,Y,V,20,cmap='jet') # trace les equipotentielles
plt.colorbar() # echelle de correspondance couleur-valeur de V
plt.axis('equal') # repere norme
plt.grid()
plt.show()

# visualisation du résultat avec imshow (attention orientation différente de la figure V[0][0] en haut à gauche)
"""plt.imshow(V,cmap='jet')
plt.colorbar() # echelle de correspondance couleur-valeur de V
plt.show()"""