#!/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)

# 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)

# 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])

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


# files

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)