#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 13 13:27:25 2017

@author: laro
"""

import numpy as np
from scipy import * 
from random import *
from pylab import *
import time
from numpy import cumsum
import matplotlib.pyplot as plt
from scipy import optimize

#antes correr ising_quf que genera las tiras de datos.. 
#aqui loadeamos, asegurense de tener en la misma carpeta los .npy o poner el path

Lx=32
Ly=32


dire="" #poner path de ser necesario(no estar en el .py en la misma carpeta que los datos)
meanE=np.load(dire+"mE.npy")
meanE2=np.load(dire+"mE2.npy")
meanM=np.load(dire+"mM.npy")
meanM2=np.load(dire+"mM.npy")

# =============================================================================
# stdE=np.load(dire+"stdE.npy")
# stdM=np.load(dire+"stdM.npy")
# =============================================================================
Tvect=np.load(dire+"Tvect.npy")


#si tienen problemas con hyper descomentar aqui con ctrl+5 y comentar el loop donde se contruyen tales mas abajo
# =============================================================================
# uposta=np.load(dire+"uposta.npy")
# cvposta=np.load(dire+"cvposta.npy")
# =============================================================================




#ploteamos los resultados
figure(2)
plot(Tvect,meanE)
suptitle("EvsT")
xlabel("T")
ylabel("mean Energy")

figure(3)
plot(Tvect,meanM)
suptitle("MvsT")
xlabel("T")
ylabel("mean Magn")




#calculamos cv y sus por sitio con desv estandar
#si hacen la cuenta les queda el CVtotal=varE/T**2
N=Lx*Ly
varE=(array(meanE2)-array(meanE)**2)*(N**2)
#el factor N**2 viene de que meanE es la media por sitio!
CV=varE/(array(Tvect))**2
cv_std=CV/N

figure(4)
plot(Tvect,cv_std)
suptitle("Cv por sitio")
xlabel("T")


varM=(array(meanE2)-array(meanE)**2)*(N**2)
CHI=varM/array(Tvect)
chi=CHI/N
figure(5)
plot(Tvect,chi)
suptitle("Chi por sitio")
xlabel("T")



#tambien podemos calcular cv con diferencias finitas (derivada numerica)
cv_num=np.gradient(meanE)/np.gradient(Tvect)
figure(6)
plot(Tvect,cv_num,label="con dif finitas")
plot(Tvect,cv_std,label="con varianza")
legend()
suptitle("Cv por sitio comparativo")
xlabel("T")





#calculemos symbolicamente u y cv de Onsager y campo medio
from sympy import symbols,sympify,sqrt,sin,cos,lambdify,tanh,coth,sinh,cosh,simplify
from sympy.plotting import plot3d
from sympy.tensor.array import derive_by_array

t= symbols('t',real=True, positive=True)
from scipy.special import ellipk
kpr=2*(tanh(2/t))**2-1
k=2*sinh(2/t)/(cosh(2/t))**2

from sympy.functions import hyper
import sympy as sy

#la solucion de onsager esta en terminos de integrales elípticas...
#no logre implementarlas en sympy... pero hay un plan b
#en el paper de Viswanathan ponen estas en terminos de funciones hipergeometricas, las cuales si salen facil con sypmy
elip=sy.pi/2*hyper((.5, .5), [1], k**2)
u0=- coth(2/t)*(1+2*kpr*elip/sy.pi)
cv=sy.simplify(sy.diff(u0, t))


#genero funcion lambda con funcion sympy "analitica"
from sympy.utilities.lambdify import lambdify
lam_u=lambdify(t,u0)
lam_cv=lambdify(t,cv)



#armemos un vector de Temperaturas mas sampleado para el cvposta y campo medio
Tvect2=linspace(.1,5,100)
#finalmente logro vectores u y cv , comentar si hay error
uposta=[]
cvposta=[]
for i in Tvect2:
    uposta.append(lam_u(i))
    cvposta.append(lam_cv(i))



#campo medio
q=4 #primeros vecinos
tc=q #notar que la Tc que predice campo medio es mas de 2 veces la Tc "posta"
l0=[]
for tt in Tvect2: #resulevo relacion de autoconsistencia iterativamente y genero un vector l0(t)
    def func(l0):
        q=4
        return np.tanh(q*l0/tt)
    l0.append(optimize.fixed_point(func,1))

#veamos que generamos
figure(45)
plot(Tvect2,l0,label="campo medio")
plot(Tvect,meanM,label="metropolis")
legend()
xlabel("T")
suptitle("magnetizacion vs T")


#generemos u y cv de campo medio
ucmedio=-.5*q*array(l0)**2
cvmedio=-q*array(l0)*array(np.gradient(l0)/np.gradient(Tvect2))




#grafiquemos todo
plt.figure(13)
plt.plot(Tvect2,uposta,label="u onsager")
plt.plot(Tvect2,ucmedio,label="u campo medio")
plt.plot(Tvect,meanE,label="u metropolis")
legend()
suptitle("Energia por sitio")
xlabel("T")


plt.figure(14)
plt.plot(Tvect2,array(cvposta),label="cv posta")
plt.plot(Tvect2,array(cvmedio),label="cv campo medio")
plt.plot(Tvect,cv_num,label="cv metropolis con derivada numerica")
plt.plot(Tvect,cv_std,label="cv metropolis con std")
legend()
suptitle("Cv por sitio")
xlabel("T")