import numpy as np
import numpy.random as npr
import math
import matplotlib.pyplot as plt
import time

start_time = time.time()
Lx = 20  # rozmiar boxa w gore/w dol
Ly = 20  # rozmiar boxa lewo/prawo
c = 0.5  # stezenie czastek w boxie
n = c * (Lx * Ly)
n = int(n)
shakes = 100  # ilosc krokow MCS
MCS = 100  # ilosc symulacji monte carlo
MCS_counter = MCS  # do wykresu, da sie pewnie lepiej ale dziala
Box = np.zeros((Lx, Ly))
class Point(object):
    def __init__(self, x, y):
        x = int(x)
        y = int(y)
        self.X = x
        self.Y = y
        self.xcounter = 0.0
        self.xcounter = float(self.xcounter)
        self.ycounter = 0.0
        self.ycounter = float(self.ycounter)

    def __repr__(self):
        return "%d, %d" % (self.X, self.Y)

    def move(self, dx, dy):

        self.X = self.X + dx
        self.Y = self.Y + dy
        if dx != 0:
            self.xcounter += dx
        if dy != 0:
            self.ycounter += dy

    def random_walk(self):
        walk = npr.randint(1, 5)

        if walk == 1:  # ruch w gore
            if self.X == 0 and filled_box[Lx - 1, self.Y] != 1:
                filled_box[self.X, self.Y] = 0
                self.X = Lx - 1
                filled_box[self.X, self.Y] = 1
                self.xcounter += -1

            elif self.X == 0 and filled_box[Lx - 1, self.Y] == 1:
                pass

            elif filled_box[self.X - 1, self.Y] != 1:
                filled_box[self.X, self.Y] = 0
                self.move(-1, 0)
                filled_box[self.X, self.Y] = 1

        if walk == 3:  # ruch w dol

            if self.X == Lx - 1 and filled_box[0, self.Y] != 1:
                filled_box[self.X, self.Y] = 0
                self.X = 0
                filled_box[self.X, self.Y] = 1
                self.xcounter += 1


            elif self.X == Lx - 1 and filled_box[0, self.Y] == 1:
                pass

            elif filled_box[self.X + 1, self.Y] != 1:
                filled_box[self.X, self.Y] = 0
                self.move(1, 0)
                filled_box[self.X, self.Y] = 1

        if walk == 2:  # ruch w lewo
            if self.Y == 0 and filled_box[self.X, Ly - 1] != 1:
                filled_box[self.X, self.Y] = 0
                self.Y = Ly - 1
                filled_box[self.X, self.Y] = 1
                self.ycounter += -1

            elif self.Y == 0 and filled_box[self.X, Ly - 1] == 1:
                pass

            elif filled_box[self.X, self.Y - 1] != 1:
                filled_box[self.X, self.Y] = 0
                self.move(0, -1)
                filled_box[self.X, self.Y] = 1

        if walk == 4:  # ruch w prawo
            if self.Y == Ly - 1 and filled_box[self.X, 0] != 1:
                filled_box[self.X, self.Y] = 0
                self.Y = 0
                filled_box[self.X, self.Y] = 1
                self.ycounter += 1

            elif self.Y == Ly - 1 and filled_box[self.X, 0] == 1:
                pass

            elif filled_box[self.X, self.Y + 1] != 1:
                filled_box[self.X, self.Y] = 0
                self.move(0, 1)
                filled_box[self.X, self.Y] = 1
def population(n):
    coords = []
    for j in range(0, Lx):
        for k in range(0, Ly):
            coords.append(Point(j, k))

    i = 0
    atoms = []
    while i < n:
        atoms.append(coords.pop(npr.randint(0, len(coords))))
        i += 1

    return atoms
def intothebox(atoms, box):
    for i in range(0, n):
        box[(atoms[i]).X, (atoms[i]).Y] = 1
    return box
def shake(pop):
    for i in pop:
        i.random_walk()
    return pop
R_vs_shakes = np.zeros((MCS, shakes))
while MCS >= 1:
    pop = []
    filled_box = []
    Box = np.zeros((Lx, Ly))

    pop = population(n)
    filled_box = intothebox(pop, Box)
    D = []

    for i in range(1, shakes + 1):  # robie tyle krokow monte carlo - tyle razy ruszam atomami

        dr = []
        pop = shake(pop)

        filled_box = intothebox(pop, Box)

        r = 0

        for j in range(0, n):  # licze dlugosc wektora do kwadratu dla kazdego atomu

            r = (math.sqrt((pop[j].xcounter ** 2) + (pop[j].ycounter ** 2))) ** 2
            dr.append(r)

        R_vs_shakes[MCS - 1, i - 1] = sum(dr) / len(
            dr)  # nparray w ktorym zapisuje sredni kwadrat wektora przemieszczenia dla kazdej symulacji po i krokach (numerowanie od 0 w nparray - dlatego -1)

    MCS -= 1
# tu licze srednie D dla (MCS) symulacji po (i) krokach - z R_vs_shakes jak wyzej
R_vs_shakes = np.sum(R_vs_shakes, axis=0) / np.size(R_vs_shakes, axis=0)
D_final = [0] * np.size(R_vs_shakes)
i = 0
for i in range(0, np.size(R_vs_shakes)):
    D_final[i] = 0.25 * R_vs_shakes[i] / (i + 1)

# rysuje srednie D dla ilus MCS
elapsed_time = time.time() - start_time
print elapsed_time

#zapisuje do pliku
output=dict(zip(range(1, len(D_final) + 1),D_final))

with open('099.csv', 'w') as f:
    for key in output.keys():
        f.write("%f \t %f\n" %(key, output[key]))

# plt.plot(range(1, len(D_final) + 1), D_final, 'ro')
# plt.xlabel('Ilosc krokow ')
# plt.ylabel('Wspolczynnik dyfuzji D')
# plt.title('Box ' + str(Lx) + 'x' + str(Ly) + ', ' + str(n) + ' czastek,' + str(MCS_counter) + ' symulacji,' + str(
#     shakes) + ' krokow MCS')
#
# plt.show()