%初始赋值

Ln=200;        %格点边长

L=zeros(Ln);   %格点矩阵

Q=120;         %总取向数

step_num=500;  %MC总步数

interval_save_jpg=20;         %图形存储间隔

interval_stastics=2;          %晶粒平均参数和相对密度统计间隔

stastics_data=zeros(step_num/interval_stastics,5);           %存储每interval_stastics次MCS后的平均晶粒尺寸和相对密度，存储格式为(MCS,grain count,average area,average diameter,relative density)

%烧结模拟过程参数赋值

T=1;    %温度参数

J1=1;   %晶界能

%初始结构的格点赋值

rand_l=randperm(Ln^2);

for i=1:Ln^2*V_pore

if rem(rand_l(i),Ln)==0

x=Ln;

y=fix(rand_l(i)/Ln);

else

x=rem(rand_l(i),Ln);

y=fix(rand_l(i)/Ln)+1;

end

L(x,y)=-1;  %标识孔洞区域

end

for i=(Ln^2*V_pore+1):Ln^2

if rem(rand_l(i),Ln)==0

x=Ln;

y=fix(rand_l(i)/Ln);

else

x=rem(rand_l(i),Ln);

y=fix(rand_l(i)/Ln)+1;

end

rand_Q=randperm(Q);

L(x,y)=rand_Q(1);  %标识晶粒区域

end

temp_L=zeros(Ln+2);    %标识边界区域标示为0，便于后续处理

temp_L(2:Ln+1,2:Ln+1)=L;

L=temp_L;

Ln=Ln+2;               %此时L边长Ln=Ln+2

s=[-1 -1

-1  0

-1  1

0 -1

0  1

1 -1

1  0

1  1];     %便于随即选取所选格点周围相邻的一个格点

%开始CAS模拟

for step=1:step_num

step   %显示MCS进程

rand_l=randperm(Ln^2);

for i=1:Ln^2    %随机选取格点

%if rem(i,1000)==0  %显示选取格点进程

%    i

%end

if rem(rand_l(i),Ln)==0   %把rand_l(i)转换成实际坐标L(x,y)

x=Ln;

y=fix(rand_l(i)/Ln);

else

x=rem(rand_l(i),Ln);

y=fix(rand_l(i)/Ln)+1;

end

if L(x,y)~=0  %如果格点不在边界区域则开始模拟

%---------------------如果选取点为晶粒格点---------------------

if L(x,y)~=-1  %如果选取点为晶粒格点

ss=s+[x y; x y; x y; x y; x y; x y; x y; x y];      %存储所选晶粒格点L(x,y)周围格点坐标

diff_grain_count=0;      %计数所选晶粒格点L(x,y)周围与之不同的晶粒格点数

same_grain_count=0;      %计数所选晶粒格点L(x,y)周围与之相同的晶粒格点数

diff_ss=zeros(8,3);      %存储所选晶粒格点L(x,y)周围与之不同的晶粒格点坐标及其取向度值

for ii=1:8

if L(ss(ii,1), ss(ii,2))~=0

if L(ss(ii,1), ss(ii,2))~=L(x,y) & L(ss(ii,1), ss(ii,2))~=-1

diff_grain_count=diff_grain_count+1;

diff_ss(diff_grain_count,1)=ss(ii,1);

diff_ss(diff_grain_count,2)=ss(ii,2);

diff_ss(diff_grain_count,3)=L(ss(ii,1),ss(ii,2));

end

if L(ss(ii,1), ss(ii,2))==L(x,y)

same_grain_count=same_grain_count+1;

end

if L(ss(ii,1), ss(ii,2))==-1

pore_count=pore_count+1;

end

end

end

total_grain_count=diff_grain_count+same_grain_count;

if diff_grain_count~=0  %如果所选晶粒格点L(x,y)周围有与之不同的晶粒格点

BG_energy=J1*diff_grain_count;    %BG_energy为格点所处晶界能

diff_ss_1=diff_ss(1:diff_grain_count,3);      %diff_ss_1存储被选格点周围取向度值与之不同的格点取向度值

diff_ss_2=unique(diff_ss_1);                  %去除diff_ss_1中的重复元素，并存储到diff_ss_2

rand_ll=randperm(length(diff_ss_2));

temp_Q=diff_ss_2(rand_ll(1));

change_BG_energy=J1*(total_grain_count-length(find(diff_ss_1==temp_Q)));

if change_BG_energy<=BG_energy

L(x,y)=temp_Q;

end

if change_BG_energy>BG_energy

set_probability=rand();

if exp(-(change_BG_energy-BG_energy)/T)>=set_probability

L(x,y)=temp_Q;

end

end

end  %对应于如果所选晶粒格点L(x,y)周围有与之不同的晶粒格点

end  %对应于如果选取点为晶粒格点

%---------------------如果选取点为晶粒格点---------------------

end  %对应于如果格点不在边界区域则开始模拟

end    %对应于随机选取格点

%========================================================后处理过程========================================================

%后处理1开始---------------每interval_save_jpg次MCS后存储图形矩阵---------------%

if rem(step,interval_save_jpg)==0

figure1=figure('visible','off','PaperPosition',[3.067 9.28 14.81 11.1],'PaperSize',[20.98 29.68]);

axes1 = axes('Layer','top','YDir','reverse','Parent',figure1);

axis(axes1,[0.5 Ln-2+0.5 0.5 Ln-2+0.5]);

image1 = image('CData',L(2:Ln-2+1,2:Ln-2+1),'CDataMapping','scaled','XData',[1 Ln],'YData',[1 Ln],'Parent',axes1);

if V_pore==0

jpg_name=strcat(num2str(step),'_','Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','T=',num2str(T),'_','J1=',num2str(J1),'.jpg');

else

jpg_name=strcat(num2str(step),'_','Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','Vpore=',num2str(V_pore),'_','T=',num2str(T),'_','J1=',num2str(J1),'_','J2=',num2str(J2),'.jpg');

end

saveas(image1,jpg_name,'jpg');

end

%后处理1结束---------------每interval_save_jpg次MCS后存储图形矩阵---------------%

%后处理2开始---------------每interval_stastics次MCS后统计平均晶粒尺寸，并存入stastics_data中---------------%

if rem(step,interval_stastics)==0

stastics_data(step/interval_stastics,1)=step;       %存储此次step

temp_L=L(2:Ln-1,2:Ln-1);      %去掉标示为0的边界

Q_exist=unique(temp_L);       %Q_exist存储L矩阵中仍存在的晶粒取向度Q

if Q_exist(1)==-1

Q_exist=Q_exist(2:length(Q_exist));

end

Q_length=length(Q_exist);     %Q_length为L矩阵中仍存在的晶粒取向度个数(不包括孔洞)

for qq=1:Q_length             %只统计具有在L矩阵中仍存在的取向度的晶粒的个数

temp_L=L;

temp_L(temp_L~=Q_exist(qq))=0;

temp_L=bwlabel(temp_L,8);

now_grainnum=max(max(temp_L));       %返回目前取向度为Q_exist(qq)的晶粒数目

stastics_data(step/interval_stastics,2)=stastics_data(step/interval_stastics,2)+now_grainnum;       %累加晶粒个数

end

stastics_data(step/interval_stastics,3)=(Ln-2)^2*(1-V_pore)/stastics_data(step/interval_stastics,2);       %统计晶粒平均面积

stastics_data(step/interval_stastics,4)=sqrt(stastics_data(step/interval_stastics,3));                     %统计晶粒平均直径

end

%后处理2结束---------------每interval_stastics次MCS后统计平均晶粒尺寸，并存入stastics_data中---------------%

%========================================================后处理过程========================================================

end

%结束MCS模拟

%写入文件

if V_pore==0

xls_name1=strcat('stastics_Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','T=',num2str(T),'_','J1=',num2str(J1),'.xls');

xls_name2=strcat('grainsize_Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','T=',num2str(T),'_','J1=',num2str(J1),'.mat');

else

xls_name1=strcat('stastics_Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','Vpore=',num2str(V_pore),'_','T=',num2str(T),'_','J1=',num2str(J1),'_','J2=',num2str(J2),'.xls');

xls_name2=strcat('grainsize_Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','Vpore=',num2str(V_pore),'_','T=',num2str(T),'_','J1=',num2str(J1),'_','J2=',num2str(J2),'.mat');

end

xlswrite(xls_name1,stastics_data);

save(xls_name2,'stastics_grainsize');