# -*- coding: utf-8 -*-
"""
Created on Mon Aug 16 11:40:33 2021

@author: eV
Simulation des vibrations d'une molécule de H2 (potentiel anharmonique)
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# =============================================================================
# Paramètres et fonctions utiles
# =============================================================================
epsilon= 4.52*(1.6*10**(-19))#en J
sigma=66*10**(-12) #en m
m= 1.67*10**(-27)#masse d'un atome d'hydrogène en kg

def derivEp(r):
    """Calcul de dEp/dr en J/m"""
    return 24*epsilon/r*((sigma/r)**6 - 2*(sigma/r)**(12))

# =============================================================================
# Graphe de l'énergie potentielle
# =============================================================================
listed=np.linspace(63e-12,120e-12,1000)
listeEp=4*epsilon*((sigma/listed)**12-(sigma/listed)**6)
plt.plot(listed*1e12,listeEp/1.6e-19)
plt.grid()
plt.xlabel("Distance H-H en pm")
plt.ylabel("Energie potentielle en eV")
plt.show()

# =============================================================================
# Simulation de la vibration 
# =============================================================================
d0=67*10**(-12) #distance de départ en m

def deriv(X,t):
    return np.array([X[1],-2/m*derivEp(X[0])])

tabt=np.linspace(0,1e-14,1000)

res=odeint(deriv,[d0,0],tabt)
plt.plot(tabt*1e15,res[:,0]*1e12)
plt.xlabel("t en fs")
plt.ylabel("r(t) en pm")
plt.grid()
plt.show()