import numpy as np
import matplotlib.pyplot as plt

## Parametres du probleme en USI
e=1.6e-19 # charge elementaire
eps0=8.85e-12 # permittivite dielectrique du vide
k=1/(4*np.pi*eps0)
    
## Maillage de l'espace
n = 100 # Nombre de points de discretisation selon chaque direction
# 1/2 longueur selon chaque direction en USI
Lx = 0.1 
Ly = 0.1
# Listes des xi et des yj 
xlist = np.linspace(-Lx,Lx,n)
ylist = np.linspace(-Ly,Ly,n)
# Matrice correspondant au maillage de l'espace 2D
X, Y = np.meshgrid(xlist,ylist)

## Caracteristiques du condensateur
qPi = 6e-14 # valeur de chaque charge ponctuelle
L_cond = 0.1 # Longueur des armatures du condensateur
e_cond = 0.03 # Distance entre les deux armatures
xP = np.linspace(0, L_cond,100 ) # discrétisation d'une plaque
distrib=np.zeros((2*len(xP),3)) # initialisation de la liste de la distribution de charges
# distribution des charges positives sur une armature
for i in range(len(xP)):
    distrib[i]=[-L_cond/2+xP[i],e_cond/2,qPi]
# distribution des charges négatives sur l'autre armature
for i in range(len(xP)):
    distrib[i+len(xP)]=[-L_cond/2+xP[i],-e_cond/2,-qPi]

## Expression du champ electrostatique et du potentiel crees par la distribution
Ex=0   # initialisation de la composante suivant (Ox) du champ E(M)
Ey=0   # initialisation de la composante suivant (Oy) du champ E(M)
V=0    # initialisation du potentiel V(M)
for i in range(len(distrib)):
    inv_r=1.0/((X-distrib[i][0])**2+(Y-distrib[i][1])**2)**(0.5)  # 1/PiM avec Pi le point où se trouve la charge et M le point où on calcule E
    Ex=Ex+k*distrib[i][2]*(X-distrib[i][0])*inv_r**3
    Ey=Ey+k*distrib[i][2]*(Y-distrib[i][1])*inv_r**3
    V=V+k*distrib[i][2]*inv_r
Norme=(Ex**2+Ey**2)**0.5   # norme de champ E(M)


## Représentations graphiques
plt.streamplot(X,Y, Ex, Ey, color=np.log(Norme), density=2, cmap='jet') # trace les lignes de champ E(M)
plt.contour(X,Y,V,20,cmap='jet') # trace les équipotentielles
plt.axis('equal') # repere norme
plt.grid()
plt.show()