[½»Á÷] ¡¾ÇóÖú¡¿°ïæÃÉÌØ¿¨ÂÞÔ­Àí

Originally posted by ³æ×Ó218 at 2010-09-10 18:20:38:
Ò»¡¢²ÄÁÏÌåϵÖÐÃÉÌØ¿¨ÂÞÄ£ÄâÖвÉÓÃÖØÒªÐÔ³éÑùµÄ±ØÒªÐÔ£¬²¢¾ßÌå½éÉÜÒ»ÖÖʵÏÖ·þ´Ó²¨¶û´ÄÂü·Ö²¼µÄÖØÒªÐÔ³éÑùµÄ·½·¨ºÍ²½Öè
¶þ¡¢Ó÷Ö×Ó¶¯Á¦Ñ§µÄ¹ÛµãºÍ·½·¨·ÖÎö¹ÌÌå±íÃæÄܵı¾ÖÊ
Èý¡¢Ó÷Ö×Ó¶¯Á¦Ñ§µÄ¹ÛµãºÍ·½·¨·ÖÎö£º¾§ ...

Originally posted by ³æ×Ó218 at 2010-09-10 18:20:38:
Ò»¡¢²ÄÁÏÌåϵÖÐÃÉÌØ¿¨ÂÞÄ£ÄâÖвÉÓÃÖØÒªÐÔ³éÑùµÄ±ØÒªÐÔ£¬²¢¾ßÌå½éÉÜÒ»ÖÖʵÏÖ·þ´Ó²¨¶û´ÄÂü·Ö²¼µÄÖØÒªÐÔ³éÑùµÄ·½·¨ºÍ²½Öè
¶þ¡¢Ó÷Ö×Ó¶¯Á¦Ñ§µÄ¹ÛµãºÍ·½·¨·ÖÎö¹ÌÌå±íÃæÄܵı¾ÖÊ
Èý¡¢Ó÷Ö×Ó¶¯Á¦Ñ§µÄ¹ÛµãºÍ·½·¨·ÖÎö£º¾§ ...

%³õʼ¸³Öµ
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');