from matplotlib.pyplot import *
import numpy


def moy_e(M,T):
    return 1.0/(numpy.exp(1.0/T)-1)-M/(numpy.exp(M/T)-1)
    
def var_e(M,T):
    return numpy.exp(1.0/T)/(numpy.exp(1.0/T)-1)**2-M**2*numpy.exp(-M/T)/(1-numpy.exp(-M/T))**2
    
def ecart_e(M,T):
    return numpy.sqrt(var_e(M,T))
    
def Cv(M,T):
    return var_e(M,T)/T**2
    
T = numpy.linspace(0.01,5,100)
figure()
M=1000
e1 = moy_e(M,T)
de1 = ecart_e(M,T)
plot(T,e1)
errorbar(T,e1,yerr=de1/2,fmt=None)
xlabel("T")
ylabel("e")
grid()
                     

cv1 = Cv(M,T)
figure()
plot(T,cv1)
xlabel("T")
ylabel("Cv")
grid()
                     

M=2
figure()
e1 = moy_e(M,T)
de1 = ecart_e(M,T)
plot(T,e1)
errorbar(T,e1,yerr=de1/2,fmt=None)
xlabel("T")
ylabel("e")
grid()
                      

cv1 = Cv(M,T)
figure()
plot(T,cv1)
xlabel("T")
ylabel("Cv")
grid()
                     

from particules import Particules

                

def courbe_T(N,M,nT,n1,n2):
    particules = Particules(N,M)
    T= numpy.arange(1,nT+1)*0.5
    E = numpy.zeros(nT)
    dE = numpy.zeros(nT)
    rejet = numpy.zeros(nT)
    for k in range(nT):
        print("T = %f"%T[k])
        (E[k],dE[k],rejet[k])=particules.iterations(n1,n2,T[k],"directe")
        particules.initialiser_sommes()
    figure()
    subplot(211)
    errorbar(T,E,yerr=dE,fmt=None)
    plot(T,E,'o')
    axis([0,numpy.max(T),0,numpy.max(E)*2])
    xlabel("T")
    ylabel("E")
    grid()
    subplot(212)
    plot(T,rejet,'o')
    xlabel('T')
    ylabel('rejet')
    grid()

                

N=10
M=10 
nT=20
n1=2000
n2=100
courbe_T(N,M,nT,n1,n2)
                

def courbe_T_metropolis(N,M,nT,n1,n2):
    particules = Particules(N,M)
    T= numpy.arange(1,nT+1)*0.5
    E = numpy.zeros(nT)
    dE = numpy.zeros(nT)
    for k in range(nT):
        print("T = %f"%T[k])
        particules.iterations(n1,n2,T[k],"metropolis")
        (E[k],dE[k],rejet)=particules.iterations(n1,n2,T[k],"metropolis")
        particules.initialiser_sommes()
    figure()
    errorbar(T,E,yerr=dE,fmt=None)
    plot(T,E,'o')
    axis([0,numpy.max(T),0,numpy.max(E)*2])
    xlabel("T")
    ylabel("E")
    grid()

                  

courbe_T_metropolis(N,M,nT,n1,n2)