#!/usr/bin/env python3

import random
from math import *

# size of the grid
L=24
# boundary thickness (must be even)
l0=4
# activity and interaction
z=3
J=5


# draw random dimers in a select area
def fill_dimers(mask,direction):
    if(direction=="h"):
        d=[1,0]
    else:
        d=[0,1]

    dimers=[]

    # keep track of which sites are occupied
    occupied=[]
    for i in range(L):
        occupied.append([])
        for j in range(L):
            occupied[i].append(0)

    for i in range(10000):
        # pick a random edge (indexed by its lower-left corner)
        e=None
        while(e==None or e[0]+d[0]>=L or e[1]+d[1]>=L or mask[e[0]][e[1]]==0 or mask[e[0]+d[0]][e[1]+d[1]]==0):
            e=[random.randint(0,L-1),random.randint(0,L-1)]
        # check whether a dimer can be added to the edge
        if(occupied[e[0]][e[1]]==0 and occupied[e[0]+d[0]][e[1]+d[1]]==0):
            # number of interactions
            interactions=0
            if(e[0]+2*d[0]<L and e[1]+2*d[1]<L and occupied[e[0]+2*d[0]][e[1]+2*d[1]]==1):
                interactions=interactions+1
            if(e[0]-d[0]>=0 and e[1]-d[1]>=0 and occupied[e[0]-d[0]][e[1]-d[1]]==1):
                interactions=interactions+1
            # probability of adding the dimer
            p=1/(1+1/z*exp(-J*interactions))
            if(p>random.random()):
                # add dimer
                dimers.append(e)
                occupied[e[0]][e[1]]=1
                occupied[e[0]+d[0]][e[1]+d[1]]=1

    return(dimers)

# find interactions
def interactions(dimers,direction):
    if(direction=="h"):
        d=[1,0]
    else:
        d=[0,1]
    
    out=[]
    for d1 in dimers:
        for d2 in dimers:
            if(d1[0]-d2[0]==2*d[0] and d1[1]-d2[1]==2*d[1]):
                out.append([d2[0]+d[0],d2[1]+d[1]])
    return(out)


# draw dimers
def draw_dimers(v_dimers,h_dimers,file_desc,color):
    for d in v_dimers:
        print("\\dimer{[color="+color+"]("+str(d[0])+","+str(d[1])+")}v", file=file_desc)
    print("", file=file_desc)
    for d in h_dimers:
        print("\\dimer{[color="+color+"]("+str(d[0])+","+str(d[1])+")}h", file=file_desc)
    print("", file=file_desc)

# draw interactions
def draw_interactions(v_interactions,h_interactions,file_desc):
    for d in v_interactions:
        print("\\interaction{("+str(d[0])+","+str(d[1])+")}v", file=file_desc)
    print("", file=file_desc)
    for d in h_interactions:
        print("\\interaction{("+str(d[0])+","+str(d[1])+")}h", file=file_desc)
    print("", file=file_desc)

# draw loops
def draw_loops(loops,file_desc,color):
    for loop in loops:
        print("\\draw[color="+color+", line width=3pt]",end="",file=file_desc)
        for e in loop:
            print("("+str(e[0])+","+str(e[1])+")--",end="",file=file_desc)
        print("cycle;",file=file_desc)

# init tikz file
def init_tikz(filename):
    file_desc=open(filename,"w")
    print("\\documentclass{standalone}\n\n\\usepackage{tikz}\n\\usepackage{dimer}\n\\usetikzlibrary{decorations.pathmorphing}\n\n\\begin{document}\n\\begin{tikzpicture}\n\n", file=file_desc)
    return(file_desc)

# close tikz file
def close_tikz(file_desc):
    print("\\end{tikzpicture}\n\\end{document}", file=file_desc)
    file_desc.close()


mu=int(L/2)

# loops
loop1=[]
for i in range(0,2*l0):
    loop1.append([mu-l0+i,mu+4.5])
for i in range(0,4):
    loop1.append([mu+l0-1+0.5,mu+4-i])
for i in range(0,l0):
    loop1.append([mu+l0+i,mu+0.5])
for i in range(0,6):
    loop1.append([mu+2*l0-1+0.5,mu-i])
for i in range(0,4*l0):
    loop1.append([mu+2*l0-1-i,mu-5.5])
for i in range(0,6):
    loop1.append([mu-2*l0-0.5,mu-5+i])
for i in range(0,l0):
    loop1.append([mu-2*l0+i,mu+0.5])
for i in range(0,4):
    loop1.append([mu-l0-0.5,mu+1+i])
# core of loop1
core_loop1=[]
for i in range(0,2*l0-4):
    core_loop1.append([mu-l0+2+i,mu+3.5])
for i in range(0,3):
    core_loop1.append([mu+l0-3+0.5,mu+3-i])
core_loop1.append([mu+l0-2,mu+0.5])
core_loop1.append([mu+l0-1,mu+0.5])
for i in range(0,l0-2):
    core_loop1.append([mu+l0+i,mu-1+0.5])
for i in range(0,4):
    core_loop1.append([mu+2*l0-3+0.5,mu-1-i])
for i in range(0,4*l0-4):
    core_loop1.append([mu+2*l0-3-i,mu-4.5])
for i in range(0,4):
    core_loop1.append([mu-2*l0+2-0.5,mu-4+i])
for i in range(0,l0-2):
    core_loop1.append([mu-2*l0+2+i,mu-1+0.5])
core_loop1.append([mu-l0,mu+0.5])
core_loop1.append([mu-l0+1,mu+0.5])
for i in range(0,3):
    core_loop1.append([mu-l0+2-0.5,mu+1+i])

loop2=[[mu-1,mu+2.5],[mu-0.5,mu+2],[mu-0.5,mu+1],[mu-1,mu+0.5],[mu-1.5,mu+1],[mu-1.5,mu+2]]
loop3=[[mu-1,mu-1.5],[mu,mu-1.5],[mu+0.5,mu-2],[mu+0.5,mu-3],[mu,mu-3.5],[mu-1,mu-3.5],[mu-1.5,mu-3],[mu-1.5,mu-2]]


loops=[loop1,loop2,loop3]


# masks
# init
v_mask=[]
h_mask=[]
for i in range(L):
    v_mask.append([])
    h_mask.append([])
    for j in range(L):
        v_mask[i].append(1)
        h_mask[i].append(0)
# draw masks
for i in range(mu-l0,mu+l0):
    for j in range(mu+1,mu+5):
        v_mask[i][j]=0
for i in range(mu-2*l0,mu+2*l0):
    for j in range(mu-5,mu+1):
        v_mask[i][j]=0
for i in range(mu-l0+2,mu+l0-2):
    for j in range(mu+1,mu+4):
        h_mask[i][j]=1
for i in range(mu-l0,mu+l0):
    h_mask[i][mu]=1
for i in range(mu-2*l0+2,mu+2*l0-2):
    for j in range(mu-4,mu):
        h_mask[i][j]=1
h_mask[mu-1][mu+1]=0
h_mask[mu-1][mu+2]=0
h_mask[mu-1][mu-3]=0
h_mask[mu-1][mu-2]=0
h_mask[mu][mu-3]=0
h_mask[mu][mu-2]=0

# random dimers in mask
v_dimers=fill_dimers(v_mask,"v")
h_dimers=fill_dimers(h_mask,"h")

# mantle dimers
h_mantle=[]
for i in range(0,l0):
    h_mantle.append([mu-l0+2*i,mu+4])
for i in range(0,4):
    h_mantle.append([mu+l0-2,mu+4-i])
    h_mantle.append([mu-l0,mu+4-i])
for i in range(0,int(l0/2)):
    h_mantle.append([mu+l0+2*i,mu])
    h_mantle.append([mu-l0-2-2*i,mu])
for i in range(0,5):
    h_mantle.append([mu+2*l0-2,mu-i])
    h_mantle.append([mu-2*l0,mu-i])
for i in range(0,2*l0):
    h_mantle.append([mu-2*l0+2*i,mu-5])

v_mantle=[]
v_mantle.append([mu-1,mu+1])
v_mantle.append([mu-1,mu-3])
v_mantle.append([mu,mu-3])

# flip horizontal mantle
flip_mantle=[]
for i in range(2,2*l0-2):
    flip_mantle.append([mu-l0+i,mu+4])
for i in range(0,2):
    flip_mantle.append([mu+l0-2,mu+3-2*i])
    flip_mantle.append([mu-l0+1,mu+3-2*i])
flip_mantle.append([mu+l0-1,mu+1])
flip_mantle.append([mu-l0,mu+1])
for i in range(1,l0-1):
    flip_mantle.append([mu+l0-1+i,mu])
    flip_mantle.append([mu-l0-i,mu])
for i in range(0,2):
    flip_mantle.append([mu+2*l0-2,mu-2*i-2])
    flip_mantle.append([mu-2*l0+1,mu-2*i-2])
for i in range(2,4*l0-2):
    flip_mantle.append([mu-2*l0+i,mu-6])

# interactions
v_interactions=interactions(v_dimers+v_mantle,"v")
h_interactions=interactions(h_dimers+h_mantle,"h")


# files

grid=init_tikz("grid.tikz.tex")
print("\\grid{"+str(L-1)+"}{"+str(L-1)+"}{(0,0)}\n", file=grid)
close_tikz(grid)

dimers=init_tikz("dimers.tikz.tex")
print("\\grid{"+str(L-1)+"}{"+str(L-1)+"}{(0,0)}\n", file=dimers)
draw_dimers(v_dimers+v_mantle,h_dimers+h_mantle,dimers,"black")
close_tikz(dimers)

interaction=init_tikz("interaction.tikz.tex")
print("\\grid{"+str(L-1)+"}{"+str(L-1)+"}{(0,0)}\n", file=interaction)
draw_interactions(v_interactions,h_interactions,interaction)
draw_dimers(v_dimers+v_mantle,h_dimers+h_mantle,interaction,"black")
close_tikz(interaction)

boundary=init_tikz("boundary.tikz.tex")
print("\\fill[color=cyan](-0.5,-0.5)--++("+str(L)+",0)--++(0,"+str(l0)+")--++("+str(-L)+",0)--cycle;", file=boundary)
print("\\fill[color=cyan](-0.5,"+str(L-l0-0.5)+")--++("+str(L)+",0)--++(0,"+str(l0)+")--++("+str(-L)+",0)--cycle;", file=boundary)
print("", file=boundary)
print("\\grid{"+str(L-1)+"}{"+str(L-1)+"}{(0,0)}\n", file=boundary)
draw_dimers(v_dimers+v_mantle,h_dimers+h_mantle,boundary,"black")
close_tikz(boundary)

contour=init_tikz("dimer_contour.tikz.tex")
print("\\grid{"+str(L-1)+"}{"+str(L-1)+"}{(0,0)}\n", file=contour)
draw_loops([loop1,loop2,loop3],contour,"blue")
draw_interactions(v_interactions,h_interactions,contour)
draw_dimers(v_dimers,h_dimers,contour,"black")
draw_dimers(v_mantle,h_mantle,contour,"black")
close_tikz(contour)