import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt

def warPocz(N):
    """generator konfiguracji spinów dla warunków początkowych"""
    spin = 2 * np.random.randint(2, size=(N,N)) - 1
    return spin

def metropolis(spin, beta):
    """Algorytm Metropolisa"""
    for i in range(N):
        for j in range (N):
            a = np.random.randint(0, N)
            b = np.random.randint(0, N)
            s = spin[a, b]
            nb = spin[(a+1)%N, b] + spin[a, (b+1)%N] + spin[(a-1)%N, b] + spin[a, (b-1)%N]
            cost = 2 * s * nb
            if cost < 0:
                s *= -1
            elif rand() < np.exp(-cost * beta):
                s *= -1
            spin[a, b] = s
    return spin

def energia(spin):
    energia = 0
    for i in range(len(spin)):
        for j in range(len(spin)):
            S = spin[i, j]
            nb = spin[(i+1)%N, j] + spin[i, (j+1)%N] + spin[(i-1)%N,j] + spin[i, (j-1)%N]
            energia += -nb * S
    return energia/4

def magnetyzacja(spin):
    mag = np.sum(spin)
    return mag

if __name__ == '__main__':
    #parametry symulacji
    nt      = 32        #temperatura
    N       = 10        #wymiar siatki 10x10
    amPowt  = 1000      #ilosc powtorzen algorytmu Metropolisa

    T = np.linspace(1.53, 3.28, nt);
    E, M, C, X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
    n1, n2 = 1.0/(amPowt*N*N), 1.0/(amPowt*amPowt*N*N)

    for t in range(nt):
        E1 = M1 = E2 = M2 = 0
        spin = warPocz(N)
        iT = 1.0/T[t]
        iT2 = iT*iT

        for i in range(amPowt):
            metropolis(spin, iT)

        for i in range(amPowt):
            metropolis(spin, iT)
            kalkEne = energia(spin)
            kalkMag = magnetyzacja(spin)
            E1 = E1 + kalkEne
            M1 = M1 + kalkMag
            E2 = E2 + kalkEne * kalkEne
            M2 = M2 + kalkMag * kalkMag

        E[t] = n1 * E1
        M[t] = n1 * M1
        C[t] = (n1 * E2 - n2 * E1 * E1) * iT2
        X[t] = (n1 * M2 - n2 * M1 * M1) * iT

# wykresy
    f = plt.figure(figsize=(18, 10));

    sp = f.add_subplot(2, 2, 1);
    plt.scatter(T, E, s=50, marker='o', color='red')
    plt.xlabel("Temperatura", fontsize=20);
    plt.ylabel("Energia", fontsize=24);
    plt.axis('tight')

    sp = f.add_subplot(2, 2, 2);
    plt.scatter(T, abs(M), s=50, marker='o', color='Blue')
    plt.xlabel("Temperatura", fontsize=20);
    plt.ylabel("Magnetyzacja", fontsize=24);
    plt.axis('tight')

    sp = f.add_subplot(2, 2, 3);
    plt.scatter(T, C, s=50, marker='o', color='purple')
    plt.xlabel("Temperatura", fontsize=20);
    plt.ylabel("Cieplo", fontsize=24);
    plt.axis('tight')

    sp = f.add_subplot(2, 2, 4);
    plt.scatter(T, X, s=50, marker='o', color='green')
    plt.xlabel("Temperatura", fontsize=20);
    plt.ylabel("Podatnosc",fontsize=24);
    plt.axis('tight')
    plt.show()