text: import numpy
  import matplotlib.pyplot as plt
  import h5py
  
  def ising(N, T, B):
  	# Initialization
  	n = 50000000
  	xs = [1, 0, -1, 0]
  	ys = [0, 1, 0, -1]
  
  	M = 0 # Magnetization
  	E = 0 # Energy
  
  	J = 0
  	k = 1
  	
  	beta = 1.0/numpy.float32(k*T)
  
  	grid = numpy.random.random_integers(0, 1, (N, N))
  	
  	for i in range(N):
  		for j in range(N):
  			if(grid[i, j] == 0):
  				grid[i, j] -= 1
  
  	# Initial Energy
  	for i in range(N):
  		for j in range(N):
  			ssum = 0
  			for x in range(4):
  				i2 = i + xs[x]
  				j2 = j + ys[x]
  				
  				if(i2 < 0):
  					i2 = N-1
  				elif(i2 >= N):
  					i2 = 0
  				if(j2 < 0):
  					j2 = N-1
  				elif(j2 >= N):
  					j2 = 0
  				ssum += grid[i2, j2]
  			E += -J * grid[i, j] * (ssum)/2 - B * grid[i, j]
  	# Initial Magnetization
  	for i in range(N):
  		for j in range(N):
  			M += grid[i, j]
  
  	# Simulation
  	for z in range(n):
  		#T = E
  		i = numpy.random.random_integers(0, N-1, 1)[0]
  		j = numpy.random.random_integers(0, N-1, 1)[0]
  
  		ssum = 0
  		
  		for x in range(4):
  			i2 = i + xs[x]
  			j2 = j + ys[x]
  			if(i2 < 0):
  				i2 = N-1
  			elif(i2 >= N):
  				i2 = 0
  			if(j2 < 0):
  				j2 = N-1
  			elif(j2 >= N):
  				j2 = 0
  			ssum += grid[i2, j2]
  
  		delta_e = 2*J*grid[i, j]*(ssum) + 2 * B * grid[i, j]
  		
  		if(delta_e <= 0):
  			grid[i, j] = -grid[i, j]
  			E += delta_e
  			M += 2*grid[i, j]
  	
  		elif(numpy.random.random(1)[0] < numpy.exp(-beta*delta_e)):#exps[delta_e]): #numpy.exp(-k*T*delta_e)):
  			grid[i, j] = -grid[i, j]
  			E += delta_e
  			M += 2*grid[i, j]
  	return (M, E)
  	
  def metropolis(B, n, T1, T2):
  	N = 50
  	T = T1
  	step = numpy.float32((T2-T1)/n)
  	print step
  	
  	print('b='+str(B)+',j=0,n=50,t='+str(T1)+'-'+str(T2)+'.h5')
  	h5f = h5py.File('N='+str(N)+',b='+str(B)+',j=0,n=50,t='+str(T1)+'-'+str(T2)+'.h5', 'w')
  
  	Ms = numpy.empty([n], dtype='f')
  	Es = numpy.empty([n], dtype='f')
  	Ts = numpy.empty([n], dtype='f')
  	
  	for i in range(n):
  		#T = numpy.random.random(1)[0]*5
  		res = ising(N, T, B)
  		Ms[i] = res[0]/numpy.float32(N*N)
  		Es[i] = res[1]/numpy.float32(N*N)
  		Ts[i] = T
  		T += step
  		print i
  
  	h5f.create_dataset('Ms', data=Ms)
  	h5f.create_dataset('Es', data=Es)
  	h5f.create_dataset('Ts', data=Ts)
  	h5f.close()
  
  
  	# h5f = h5py.File('data.h5','r')
  	# b = h5f['dataset_1'][:]
  	# h5f.close()
  
  	# plt.figure()
  	# plt.ylabel('Magnetization per site')
  	# plt.xlabel('Temperature')
  	# plt.scatter(Ts, Ms)
  	# plt.savefig('b='+str(5.0)+',j=1,n=50,t=2.0.png')
  	# plt.figure()
  	# plt.ylabel('Energy per site')
  	# plt.xlabel('Temperature')
  	# plt.scatter(Ts, Es)
  	# plt.figure()
  	# plt.ylabel('Magnetization per site')
  	# plt.xlabel('Energy per site')
  	# plt.scatter(Es, Ms)
  	#plt.axis([0, 5, -3.5, 3.5])	
  	#plt.show()
  
  
  
  metropolis(0.3, 20, 0.1, 2.0)