#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 18 10:39:13 2022

@author: administrateur
"""

import numpy as np
import matplotlib.pyplot as plt

plt.close('all')

### Constantes physiques
eps0 = 8.85e-12  # en F.m-1

### Caractéristiques de la molécule
a = 127.4   # en pm
q = 2.9e-20 # en coulomb

x_Cl = 0
y_Cl = -a/2

x_H = 0
y_H = a/2

### Caractéristiques de la grille
h = 150     # hauteur de la grille en pm
Nx = 101    # nbre de pts
Ny = 101
dx = h/Nx   # pas d'espace 
dy = h/Ny


### Construction de la grille
xx = np.linspace(-h/2 , h/2, num=Nx) # grille symétrique
yy = np.linspace(-h/2 , h/2, num=Ny)
x, y = np.meshgrid(xx, yy)


### Calcul du potentiel :
V = np.empty((Nx,Ny))
for i in range(len(xx)):
    for j in range(len(yy)):
        r_H = np.sqrt( (x[i,j]-x_H)**2 + (y[i,j]-y_H)**2 )
        r_Cl = np.sqrt( (x[i,j]-x_Cl)**2 + (y[i,j]-y_Cl)**2 )
        if r_H == 0 or r_Cl == 0:
            V[i,j] = 0
        else:
            V[i,j] = q/(4*np.pi*eps0) * (1/r_H - 1/r_Cl)
        

### Calcul du champ
Ey, Ex = np.gradient(-V)
### np.gradient renvoie une **liste** dont les éléments sont des np.array,
### le signe - est donc à mettre à l'intérieur de la fonction


### Tracés :

# Ce qui est important ...
plt.figure()
plt.contour(x,y,V, 500)   # 500 équipotentielles
plt.streamplot(x,y,Ex,Ey, color='k', linewidth=.5)

# ... et de quoi décorer
plt.plot([0],[-a/2], "yo", markersize = 20, markeredgecolor='k')
plt.text(0,-a/2,'Cl', ha='center', va='center')
plt.plot([0],[a/2], "yo", markersize = 20, markeredgecolor='k')
plt.text(0,a/2,'H', ha='center', va='center')
plt.axis('scaled')
plt.xlim(np.min(xx),np.max(xx))
plt.ylim(np.min(yy),np.max(yy))


# plt.savefig('./12_potentiel_fig-simu-HCl.pdf', bbox_inches='tight')