from matplotlib.pyplot import *
import math
import numpy

def force(r):
    return 48.0*(math.pow(r,-13)-0.5*math.pow(r,-7))
r = numpy.arange(0.8,4.0,0.01)
f = r.copy()
for i in range(r.size):
    f[i] = force(r[i])
figure(figsize=(8,6))
plot(r,f)
xlabel("r/d")
ylabel("force")
axis([0,4,-3,10])
grid()
                  

from lennardJones2d import *
import pygame
from matplotlib.pyplot import *
import time
            

Nx = 20
densite = 0.3
rc = 2.5
h = 0.01
sys = Systeme(Nx,densite,rc)
sys.initialiser(1.0)
pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/sys.L
clock = pygame.time.Clock()
done = False
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
iter = 500
t = time.clock()
while not done and iter > 0:
    iter -= 1
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            done = True
    sys.dessiner_disques(screen,echelle,(255,0,0))
    sys.integration(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    pygame.display.flip()
print(time.clock()-t)
pygame.image.save(screen,"../../../../figures/sciphys/dynmol/lennardjones2d/disques1.png")
pygame.quit()

           

figure(figsize=(10,4))
plot(energie)
axis([0,energie.size,-2.0,2.0])
title("Energie")
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Energie cinetique")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           

Nx = 20
densite = 0.3
rc = 2.5
deltaR=0.5
h = 0.01
sys = Systeme(Nx,densite,rc,deltaR)
sys.initialiser(1.0)
pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/sys.L
clock = pygame.time.Clock()
done = False
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
iter = 500
t = time.clock()
while not done and iter > 0:
    iter -= 1
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            done = True
    sys.dessiner_disques(screen,echelle,(255,0,0))
    sys.integration_avec_liste_voisins(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    pygame.display.flip()
print(time.clock()-t)
pygame.image.save(screen,"../../../../figures/sciphys/dynmol/lennardjones2d/disques2.png")
pygame.quit()

           

figure(figsize=(10,4))
plot(energie)
axis([0,energie.size,-2.0,2.0])
title("Energie")
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           

Nx = 40
densite = 0.3
rc = 2.5
deltaR=0.5
h = 0.01
sys = Systeme(Nx,densite,rc,deltaR)
sys.initialiser(1.0)
pygame.init()
taille = 500
screen = pygame.display.set_mode([taille,taille])
echelle = taille*1.0/sys.L
clock = pygame.time.Clock()
done = False
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
iter = 500
t = time.clock()
while not done and iter > 0:
    iter -= 1
    clock.tick(30)
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            done = True
    sys.dessiner_disques(screen,echelle,(255,0,0))
    sys.integration_avec_liste_voisins(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    pygame.display.flip()
print(time.clock()-t)
pygame.image.save(screen,"../../../../figures/sciphys/dynmol/lennardjones2d/disques3.png")
pygame.quit()

           

figure(figsize=(10,4))
plot(energie)
title("Energie")
axis([0,energie.size,-2.0,2.0])
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,0,pression.max()])
title("pression")
           

Nx = 40
densite = 0.3
rc = 2.5
deltaR=0.5
h = 0.01
T = 0.5
vinit = math.sqrt(2*T)
pression = numpy.zeros(0)
ec = numpy.zeros(0)
energie = numpy.zeros(0)
sys = Systeme(Nx,densite,rc,deltaR)
sys.initialiser(vinit)
sys.integration_avec_liste_voisins(h,1000)
sys.init_moyennes()
iter = 1000
while iter>0:
    iter -= 1
    sys.integration_avec_liste_voisins(h,10)
    sys.calculer_cinetique()
    pression = numpy.append(pression,sys.pression)
    ec = numpy.append(ec,sys.Ecinetique)
    energie = numpy.append(energie,sys.energie)
    if iter%100==0:
        sys.ajuster_vitesses(T)
            

figure(figsize=(10,4))
plot(energie)
title("Energie")
axis([0,energie.size,-2.0,2.0])
           

figure(figsize=(10,4))
plot(ec)
axis([0,ec.size,0,ec.max()])
title("Temperature")
           

figure(figsize=(10,4))
plot(pression)
axis([0,pression.size,-1,1])
title("pression")