import numpy as np
import matplotlib.pyplot as plt

n = 100 # n+1 points de discretisation sur l'axe (Ox)
m = 100 # m+1 points de discretisation sur l'axe (Oy)
# Données numériques : :

Delta_x = 0.001 # pas de discrétisation en m (pour mémoire)
L = 40. # largeur des armatures en unité de Delta_x
a = 20. # distance entre les armatures en unité de Delta_x
e = 6. # épaisseur des armatures en unités de Delta_x
U = 10. # potentiel de référence en V

# Paramètres du calcul :
eps = 1e-4*U # choix de la précision
nbMaxIter = 2000 # nombre maximum d'itérations

# Initialisation des tableaux des potentiels :
V0 = np.zeros((n+1,m+1)) # création d'un tableau nul
V1 = np.zeros((n+1,m+1)) # création d'un tableau nul

# Conditions aux limites :
def condLim(V) : # cette fonction modifie le tableau V pour qu'il vérifie les conditions aux limites :
# a) nullité de V sur la frontière du domaine de calcul :
    for i in range(n+1) :
        V[i,0] = 0
        V[i,m] = 0
    for j in range(m+1) :
        V[0,j] = 0
        V[n,j] = 0
    # b) valeurs du potentiel imposées sur les rubans conducteurs :
    imin, imax = n//2-int(L)//2, n//2+int(L)//2 #valeurs extrêmes de i pour les rubans
    jminS, jmaxS = m//2 + int(a-e)//2, m//2 + int(a+e)//2 # valeurs extrêmes de j pour le ruban supérieur
    jminI, jmaxI = m//2 + int(-a-e)//2, m//2 + int(-a+e)//2 # valeurs extrêmes de j pour le ruban inférieur

    for i in range(imin, imax+1) :
        for j in range(jminS,jmaxS+1) :
            V[i,j] = U
        for j in range(jminI,jmaxI+1) :
            V[i,j] = -U

# Methode de relaxation :
def distance(V0, V1) : # calcul de l'écart entre les deux tableaux de potentiels
    d = 0
    for i in range(n+1) :
        for j in range(m+1) :
            d = max(d, abs(V1[i,j]- V0[i,j]))
    return d


d = np.inf # réel flottant supérieur à tous les autres
iter = 0
# Écrire ci-dessous la boucle correspondant à la méthode de relaxation


# Représentation graphique :
x = np.array([i*Delta_x for i in range(n+1)]) # tableau des abscisses
y = np.array([j*Delta_x for j in range(m+1)]) # tableau des ordonnées
valeurs = [-U + k*2*U/10 for k in range(1,11)] # liste des valeurs des équipotentielles à tracer
plt.contour(x,y,np.transpose(V1),valeurs,colors = 'gray') # tracé des équipotentielles (il faut donner à cette fonction la transposée de V1, cf. https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.contour.html)

def traceRectangle(xmin,xmax,ymin,ymax) : # fonction qui trace le rectangle de sommets opposés (xmin, ymin) et (xmax, ymax)
    plt.plot((xmin,xmin,xmax,xmax, xmin),(ymin, ymax,ymax, ymin,ymin
),color='black')

traceRectangle(Delta_x*(n/2-L/2), Delta_x*(n/2+L/2), Delta_x*(m/2+a
/2-e/2), Delta_x*(m/2+a/2+e/2)) # dessin de l'armature supérieure
traceRectangle(Delta_x*(n/2-L/2), Delta_x*(n/2+L/2), Delta_x*(m/2-a
/2-e/2), Delta_x*(m/2-a/2+e/2)) # dessin de l'armature inférieure
plt.axis('equal') # pour avoir la même échelle sur les deux axes
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()