import numpy
import random
import math
from matplotlib.pyplot import *
import ising1D
                

N=1000
ising = ising1D.Ising1D(N)
ising.constantes(0.0,1.0)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()
                

ising.constantes(0.0,10.0)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()
                

ising.constantes(0.0,0.1)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid() 
                

ising.constantes(5.0,10.0)
(m,M,dM)=ising.boucle(100)
figure(figsize=(10,6))
plot(m)
xlabel('iteration')
ylabel('M')
axis([0,100,-N,N])
grid()
                

bmax = 20.0
B = numpy.arange(0.0,bmax,bmax/10)
n = B.size
listM = numpy.zeros(n)
listDM = numpy.zeros(n)
T = 10.0
N=1000
ising = ising1D.Ising1D(N)
ising.constantes(0.0,T)
ising.boucle(50)
for k in range(n):
    ising.constantes(B[k],T)
    ising.boucle(10)
    (m,M,dM)=ising.boucle(100)
    listM[k] = M
    listDM[k] = dM
figure(figsize=(8,6))
errorbar(B,listM,yerr=listDM,fmt=None)
axis([0,bmax,0,N])
ylabel("M")
xlabel("B")
title("T=10")
grid()
def moment(B,T):
    sh = math.sinh(B/(1*T))
    return N*sh/math.sqrt(sh*sh+math.exp(-4.0/(T)))
B = numpy.arange(0.0,bmax,bmax/100)
M = numpy.zeros(B.size)
for k in range(B.size):
    M[k] = moment(B[k],T)
plot(B,M)
                

bmax = 20.0
B = numpy.arange(0.0,bmax,bmax/10)
n = B.size
listM = numpy.zeros(n)
listDM = numpy.zeros(n)
T = 5.0
N=1000
ising = ising1D.Ising1D(N)
ising.constantes(0.0,T)
ising.boucle(50)
for k in range(n):
    ising.constantes(B[k],T)
    ising.boucle(10)
    (m,M,dM)=ising.boucle(100)
    listM[k] = M
    listDM[k] = dM
figure(figsize=(8,6))
errorbar(B,listM,yerr=listDM,fmt=None)
axis([0,bmax,0,N])
ylabel("M")
xlabel("B")
title("T=5")
grid()
B = numpy.arange(0.0,bmax,bmax/100)
M = numpy.zeros(B.size)
for k in range(B.size):
    M[k] = moment(B[k],T)
plot(B,M)