# -*- coding: utf-8 -*-
# Simulação de Monte Carlo para uma coleção de N sistemas de dois níveis não interagentes
# A energia é medida em unidades de épsilon, a diferença entre os dois níveis
from __future__ import division, print_function
from visual import *
from visual.graph import *
import numpy as np
    
# Inicializa os parâmetros
N = 100                    # Número de sistemas
passos = 10                # Número total de passos de passos de Monte Carlo por sistema
Tmin = 0.1                 # Temperatura mínima
Tmax = 2.                  # Temperatura máxima
dT = 0.1                   # Incremento de temperatura
# Cria o esqueleto da visualização
g_energia = gdisplay(title="Energia vs. T",ymax=0.5)           # Tela que registrará a energia
g_calesp = gdisplay(title="Calor específico vs. T",ymax=0.6)   # Tela que registrará o calor específico
energia_exata = gcurve(color=color.blue,gdisplay=g_energia)    # Curva da energia exata em azul
energia_simul = gdots(color=color.cyan,gdisplay=g_energia)     # Pontos da energia simulada em ciano
calesp_exato = gcurve(color=color.red,gdisplay=g_calesp)       # Curva do calor específico exato em vermelho
calesp_simul = gdots(color=color.magenta,gdisplay=g_calesp)    # Pontos do calor específico simulado em magenta

# Traça as curvas exatas da energia e do calor específico (por sistema) em função da temperatura
T = Tmin
while T<Tmax:
	boltz = exp(-1./T) # Fator de Boltzmann
	uboltz = 1.+boltz
	E = boltz/uboltz
	C = boltz*((1./T/uboltz)**2)
	energia_exata.plot(pos=(T,E))
	calesp_exato.plot(pos=(T,C))
	T = T+dT/10.       # Incrementa a temperatura

# Inicializando os vetores e matrizes; definindo as condições iniciais
iTmax = 1+int((Tmax-Tmin)/dT)
energia = [[0. for iT in arange(iTmax+1)] for i in arange(N+1)] # Estado de cada sistema (inicial = 0) em cada temperatura
# Vetores da simulação
energia_acum = [0. for iT in arange(iTmax+1)]                   # Acumulador da energia em cada temperatura
energia2_acum = [0. for iT in arange(iTmax+1)]                  # Acumulador do quadradado da energia em cada temperatura
energia_media = [0. for iT in arange(iTmax+1)]                  # Energia média, por sistema, em cada temperatura
calesp_medio = [0. for iT in arange(iTmax+1)]                   # Calor específico médio em cada temperatura
# Vetores do cálculo exato
energia_ex = [0. for iT in arange(iTmax+1)]                     # Energia interna exata, por sistema, em cada temperatura
calesp_ex = [0. for iT in arange(iTmax+1)]                      # Calor específico exato em cada temperatura

# Registrando os valores exatos nas mesmas temperaturas da simulação,
# para determinação posterior do erro
iT = 0
while iT<=iTmax:
	T = Tmin+iT*dT
	boltz = exp(-1./T)
	uboltz = 1.+boltz
	energia_ex[iT] = boltz/uboltz
	calesp_ex[iT] = boltz*((1./T/uboltz)**2)
	iT = iT+1

# Início da simulação; laço sobre os passos de MC
m=1
while m<=passos:
	erro_energia = 0.
	erro_calesp = 0.
    # Laço sobre as temperaturas
	iT = 0
	while iT<=iTmax:
		rate(1000)
		T = Tmin+iT*dT
		# Varrendo os sistemas (sequencialmente)
		i=1
		while i<=N:
			energia[i][iT] = 1.-energia[i][iT]     # Invertemos o estado do sistema
			delta_E = 2.*energia[i][iT]-1.         # Variação da energia na inversão
			boltz = exp(-delta_E/T)                # Fator de Boltzmann da inversão
			p = random.random()                    # Sorteamos um número aleatório entre 0 e 1
			if p > boltz:                          # Se a inversão é rejeitada,
				energia[i][iT] = 1.-energia[i][iT] # retornamos a energia do sistema ao valor inicial
			i = i+1                                # Passamos para o próximo sistema  
		# Após um passo de Monte Carlo por sistema em dada temperatura, atualizamos as medidas
		soma_energia = 0.
		i=1
		while i<=N:  # Calcula a energia total em uma certa configuração
			soma_energia = soma_energia + energia[i][iT]
			i = i+1
		energia_acum[iT] = energia_acum[iT]+soma_energia                        # Acumula a energia
		energia2_acum[iT] = energia2_acum[iT]+soma_energia**2                   # Acumula o quadrado da energia
		energia_media[iT] = energia_acum[iT]/(N*m)                              # Atualiza a energia média
		calesp_medio[iT] = (energia2_acum[iT]/m-(energia_acum[iT]/m)**2)/N/T**2 # Atualiza o calor específico médio
		# Registrando as médias nos gráficos
		energia_simul.plot(pos=(T,energia_media[iT]))
		calesp_simul.plot(pos=(T,calesp_medio[iT]))
		# Atualizando os erros e passando para a próxima temperatura
		erro_energia = erro_energia + (energia_media[iT]-energia_ex[iT])**2
		erro_calesp = erro_calesp + (calesp_medio[iT]-calesp_ex[iT])**2
		iT = iT + 1
    # "Apagando" os dados da simulação para mostrar as médias após o próximo passo
	if m<passos:
		iT = 0
		while iT<=iTmax:
			T = Tmin+iT*dT
			energia_simul.plot(pos=(T,energia_media[iT]),color=color.black)
			calesp_simul.plot(pos=(T,calesp_medio[iT]),color=color.black)
			iT = iT+1
	erro_energia = sqrt(erro_energia)/(iTmax+1)
	erro_calesp = sqrt(erro_calesp)/(iTmax+1)
	print ('Passo = ',m,' Erro na energia = ', erro_energia, ' Erro no calor esp. = ',erro_calesp)
	m = m+1