# initialiser une matrice de spins de taille h x w;
# - les elements sont +1/-1 avec une probabilite 1/2;
# - les elements sur le bord sont initialise a bc
#   ('bc' pour 'boundary conditions', doit etre 0, +1, -1)

init.ising <- function(h, w, bc=0) {
	spins <- sign(runif(h*w, -1, 1))
	M <- matrix(spins, h, w)
	M[h,] <- M[1,] <- bc
	M[,w] <- M[,1] <- bc
	return(M)
}

# visualiser une matrice de spins pour le modele d'Ising
# les spins positifs sont representes par le rouge,
# negatifs par le bleu, zeros par le gris.
# 
# parametres: 
# M   -- la matrice de spins a visualiser,
# pch -- type de point (voir "help(points)", "example(points)" 
#        pour plus d'information)

plot.ising <- function(M, pch=15) {
	w <- dim(M)[1]	
	h <- dim(M)[2]

	x <- matrix( rep(1:w, h), w, h )
	y <- t(matrix( rep(1:h, w), h, w ))

	plot(c(0,w+1,w+1,0,0), c(0, 0, h+1, h+1,0), type='l')
	points( x[M>0], y[M>0], col='red', pch=pch)
	points( x[M<0], y[M<0], col='blue', pch=pch)
	points( x[M==0], y[M==0], col='grey', pch=pch)

}


# exemple:
M <- init.ising(50, 50, 0)
plot.ising(M)