# TRACÉ D'UNE CARTE DE CHAMP ET DE POTENTIEL ELECTROSTATIQUE

import numpy as np
import matplotlib.pyplot as plt

# définition de la zone de calcul x, y
x=np.linspace(-10,10,500)
y=np.linspace(-10,10,500)
dx=20/500
# Définition des cherges et de leurs coordonnées
# on prend les charges en unités de 1/(4.Pi.epsilon0) pour éviter de trainer
# cette constante

q=[1,-3,2]
xq=[-5,5,2]
yq=[0,0,3]

# création de Ex, Ey et V nuls à la base.
Ex=np.zeros((500,500))
Ey=np.zeros((500,500))
V=np.zeros((500,500))

X,Y=np.meshgrid(x,y)                 # créé deux tableaux X,Y où on retrouve
                                     # les abcisses et ordonnées x,y et
                                     # permet un calcul vectoriel sur la grille
                                     # X, Y (cf fiche meshgrid)

# calcul des composantes de E et du potentiel V en chaque points
# pour chaque charge (il y en a len(q) le nombre d'élément du tableau q)

for i in range(len(q)):

    r=((X-xq[i])**2+(Y-yq[i])**2)**0.5 # calcul de la distance à la charge
                                       # pour chaque valeur de X et Y

    costheta=(X-xq[i])/r                     # calcul du cosinus pour la projection
    sintheta=(Y-yq[i])/r                     # et du sinus

    Ex = Ex+q[i]/r**2*costheta       # on superspose les Ex
    Ey = Ey+q[i]/r**2*sintheta       # les Ey

    V = V+q[i]/r                     # les potentiels V


# Une autre solution : on retire les calculs directs de Ex et Ey
# et on fait le calcul de -grad(V) avec une routine "maison"
#
# for i in range(1,499):
#     for j in range(1,499):
#         Ey[i,j]=-(V[i+1,j]-V[i-1,j])/(2*dx)
#         Ex[i,j]=-(V[i,j+1]-V[i,j-1])/(2*dx)

#-----------------------------------------------------------------------------
# Cette partie évite la divergence de V au voisinage des charges
# la routine "sature" V à une valeur arbitraire pour que l'affichage
# proche des charges soit visualisable.

for i in range(500):
    for j in range(500):
        if V[i,j]<-1:
            V[i,j]=-1
        if V[i,j]>1:
            V[i,j]=1
#-----------------------------------------------------------------------------

# tracé des lignes de champs : streamplot et du potentiel : contourf
# la fonction contour tracerait les équipotentielles, contourf trace en
# dégradé le potentiel V
# steamplot trace le champs de vecteur E

plt.figure(figsize=(10,8))
plt.streamplot(X,Y,Ex,Ey,density=1.5,color='black')
cf=plt.contourf(X,Y,V,100,cmap='jet')
plt.colorbar(cf)                       # échelle de couleur sur le coté
plt.show()