import numpy as np

# Parameters
J = 1
k_B = 1
T = 2.269  # Critical temperature
K = J / (k_B * T)

# Transfer matrix for a 2x2 Ising model (simplified example)
def transfer_matrix(K):
    T = np.zeros((4, 4))
    configs = [(-1, -1), (-1, 1), (1, -1), (1, 1)]
    for i, (s1, s2) in enumerate(configs):
        for j, (s3, s4) in enumerate(configs):
            E = K * (s1*s2 + s2*s3 + s3*s4 + s4*s1)
            T[i, j] = np.exp(E)
    return T

T_matrix = transfer_matrix(K)
eigenvalues, _ = np.linalg.eig(T_matrix)
lambda_max = np.max(eigenvalues)

# Free energy per spin
f = -k_B * T * np.log(lambda_max)
print(f"Free energy per spin: {f:.4f}")

# Output the transfer matrix and its largest eigenvalue
print("Transfer matrix:\n", T_matrix)
print("Largest eigenvalue:", lambda_max)

import numpy as np
import random
import math

def MC_step(config, beta):
    '''Monte Carlo move using Metropolis algorithm '''
    L = len(config)
    for _ in range(L * L):  # Corrected the loop
        a = np.random.randint(0, L)
        b = np.random.randint(0, L)
        sigma = config[a, b]
        neighbors = (config[(a + 1) % L, b] + config[a, (b + 1) % L] +
                     config[(a - 1) % L, b] + config[a, (b - 1) % L])
        del_E = 2 * sigma * neighbors
        if del_E < 0 or random.uniform(0, 1) < np.exp(-del_E * beta):
            config[a, b] = -sigma
    return config


def E_dimensionless(config, L):
    '''Calculate the energy of the configuration'''
    energy = 0
    for i in range(L):
        for j in range(L):
            S = config[i, j]
            neighbors = (config[(i + 1) % L, j] + config[i, (j + 1) % L] +
                         config[(i - 1) % L, j] + config[i, (j - 1) % L])
            energy += -neighbors * S
    return energy / 4  # To compensate for overcounting

def magnetization(config):
    '''Calculate the magnetization of the configuration'''
    return np.sum(config)


def calcul_energy_mag_C_X(config, L, eqSteps, err_runs):
    print('finished')
    
    nt = 100  # number of temperature points
    mcSteps = 1000
    T_c = 2 / math.log(1 + math.sqrt(2))
    
    T = np.linspace(1., 7., nt)
    E, M, C, X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
    C_theoric, M_theoric = np.zeros(nt), np.zeros(nt)
    delta_E, delta_M, delta_C, delta_X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
    
    n1 = 1.0 / (mcSteps * L * L)
    n2 = 1.0 / (mcSteps * mcSteps * L * L)
    
    Energies, Magnetizations, SpecificHeats, Susceptibilities = [], [], [], []
    delEnergies, delMagnetizations, delSpecificHeats, delSusceptibilities = [], [], [], []
    
    for t in range(nt):
        beta = 1. / T[t]
        
        # Equilibrate the system
        for _ in range(eqSteps):
            MC_step(config, beta)
        
        Ez, Cz, Mz, Xz = [], [], [], []
        
        for _ in range(err_runs):
            E, E_squared, M, M_squared = 0, 0, 0, 0
            
            for _ in range(mcSteps):
                MC_step(config, beta)
                energy = E_dimensionless(config, L)
                mag = abs(magnetization(config))
                
                E += energy
                E_squared += energy ** 2
                M += mag
                M_squared += mag ** 2
            
            E_mean = E / mcSteps
            E_squared_mean = E_squared / mcSteps
            M_mean = M / mcSteps
            M_squared_mean = M_squared / mcSteps
            
            Energy = E_mean / (L ** 2)
            SpecificHeat = beta ** 2 * (E_squared_mean - E_mean ** 2) / (L ** 2)
            Magnetization = M_mean / (L ** 2)
            Susceptibility = beta * (M_squared_mean - M_mean ** 2) / (L ** 2)
            
            Ez.append(Energy)
            Cz.append(SpecificHeat)
            Mz.append(Magnetization)
            Xz.append(Susceptibility)
        
        Energies.append(np.mean(Ez))
        delEnergies.append(np.std(Ez))
        Magnetizations.append(np.mean(Mz))
        delMagnetizations.append(np.std(Mz))
        SpecificHeats.append(np.mean(Cz))
        delSpecificHeats.append(np.std(Cz))
        Susceptibilities.append(np.mean(Xz))
        delSusceptibilities.append(np.std(Xz))
        
        if T[t] < T_c:
            M_theoric[t] = pow(1 - pow(np.sinh(2 * beta), -4), 1 / 8)
            C_theoric[t] = (2.0 / np.pi) * (math.log(1 + math.sqrt(2)) ** 2) * (-math.log(1 - T[t] / T_c) + math.log(1.0 / math.log(1 + math.sqrt(2))) - (1 + np.pi / 4))
        else:
            C_theoric[t] = 0
    
    return (T, Energies, Magnetizations, SpecificHeats, Susceptibilities,
            delEnergies, delMagnetizations, M_theoric, C_theoric, delSpecificHeats, delSusceptibilities)


# Parameters
L = 10  # Size of the lattice
eqSteps = 1000  # Number of steps to reach equilibrium
err_runs = 10  # Number of error runs

# Initial configuration (random spins)
config = 2 * np.random.randint(2, size=(L, L)) - 1

# Perform calculations
results = calcul_energy_mag_C_X(config, L, eqSteps, err_runs)

# Unpack results
(T, Energies, Magnetizations, SpecificHeats, Susceptibilities, 
 delEnergies, delMagnetizations, M_theoric, C_theoric, 
 delSpecificHeats, delSusceptibilities) = results

# Plot results
import matplotlib.pyplot as plt

plt.figure()
plt.errorbar(T, Energies, yerr=delEnergies, label='Energy')
plt.xlabel('Temperature')
plt.ylabel('Energy')
plt.legend()
plt.show()

plt.figure()
plt.errorbar(T, Magnetizations, yerr=delMagnetizations, label='Magnetization')
plt.xlabel('Temperature')
plt.ylabel('Magnetization')
plt.legend()
plt.show()

plt.figure()
plt.errorbar(T, SpecificHeats, yerr=delSpecificHeats, label='Specific Heat')
plt.plot(T, C_theoric, label='Theoretical Specific Heat')
plt.xlabel('Temperature')
plt.ylabel('Specific Heat')
plt.legend()
plt.show()

plt.figure()
plt.errorbar(T, Susceptibilities, yerr=delSusceptibilities, label='Susceptibility')
plt.xlabel('Temperature')
plt.ylabel('Susceptibility')
plt.legend()
plt.show()
