Monte-Carlo-Simulations/monte_carlo.m at master · nkc-137/Monte-Carlo-Simulations · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
%% To discretize voronoi diagram generated by the function
a=rand(1,300);
b=rand(1,300);
Nx = 100;
Ny = 100;
grain_old = zeros(Nx,Ny);
x = linspace(0,1,Nx);
y = linspace(0,1,Ny);
DT = delaunayTriangulation(a',b');
[v,R] = voronoiDiagram(DT);
[X,Y] = meshgrid(x,y);
G=randi([1 36],30);
for k = 1:length(R)
    if all(R{k}~=1)
        [in,on] = inpolygon(X,Y,v(R{k},1),v(R{k},2));
        for i = 1:Nx
            for j = 1:Ny
                if(in(i,j))
                    grain_old(i,j) = G(k);
                end
            end
        end
    end
end
pcolor(grain_old)
view(2);


% Initialization
size_x = 100;
size_y = 100;
grains = 30;
p = zeros(size_x,size_y);
initial_state = int16(abs(rand(size_x,size_y)*(grains-1))+1);
print_time = 1000; % Number of time steps after which the image should be printed
disp_time = 200; % Number of time steps after which the image should be shown

% p(:,:) = initial_state;
p = grain_old;
k = 1.38*10^(-23);
T = 300;
J = k*T;
N = 3000;

t = 0;
m = 1;
n = 1;

% Contains path to save the images of simulations generated
folder_name = strcat(pwd,'/images_output');

% Time loop
while (t<=10000)
    H1 = 0;
    H2 = 0;
    delta_H = 0;
    l = 0;
    while (l<N)
        p_x = int16(abs(rand(1)*(size_x-3))+2);
        p_y = int16(abs(rand(1)*(size_y-3))+2);
        H1  = H1 + J/2*(8 - kdelta(p(p_x,p_y),p(p_x+1,p_y))-kdelta(p(p_x,p_y),p(p_x-1,p_y))-kdelta(p(p_x,p_y),p(p_x,p_y+1))-kdelta(p(p_x,p_y),p(p_x,p_y-1))-kdelta(p(p_x,p_y),p(p_x+1,p_y-1))-kdelta(p(p_x,p_y),p(p_x-1,p_y-1))-kdelta(p(p_x,p_y),p(p_x-1,p_y+1))-kdelta(p(p_x,p_y),p(p_x+1,p_y+1)));

    g = int16(rand(1)*7 + 1); % To randomly select a grid around the current grid
    if (g == 1)
        m = -1;
        n = -1;
    elseif (g==2)
        m = -1;
        n = 0;
        elseif (g==3)
        m = -1;
        n = 1;
        elseif (g==4)
        m = 0;
        n = -1;
        elseif (g==5)
        m = 0;
        n = 1;
        elseif (g==6)
        m = 1;
        n = -1;
        elseif (g==7)
        m = 1;
        n = 0;
        elseif (g==8)
        m = 1;
        n = 1;
    end
    temp = p(p_x,p_y);
    p(p_x,p_y) = p(p_x+m,p_y+n);
    H2  = H2 + J/2*(8 - kdelta(p(p_x,p_y),p(p_x+1,p_y))-kdelta(p(p_x,p_y),p(p_x-1,p_y))-kdelta(p(p_x,p_y),p(p_x,p_y+1))-kdelta(p(p_x,p_y),p(p_x,p_y-1))-kdelta(p(p_x,p_y),p(p_x+1,p_y-1))-kdelta(p(p_x,p_y),p(p_x-1,p_y-1))-kdelta(p(p_x,p_y),p(p_x-1,p_y+1))-kdelta(p(p_x,p_y),p(p_x+1,p_y+1)));
    delta_H = H2-H1;

    prob = exp(-delta_H/(J));
    z = rand(1);
    if z>prob
        p(p_x,p_y) = temp;
    else
        p(p_x,p_y) = p(p_x+m,p_y+n);
    end
    l = l+1;
    end

    if (mod(t,disp_time)==0)
    pcolor(p);
    shading flat;
    view(2);
    pause(0);
    end
    if (mod(t,print_time)==0)
        disp(t);
        file_name = sprintf('image_%d',t);
        full_file_name = fullfile(folder_name,file_name);
        print(full_file_name,'-dpng');
    end
    t = t + 1;
end