%計(jì)算重心及能量
    D=zeros(222);
    E=zeros(222);
    E(:)=100;
    A1=zeros(222);
    A1(:)=100;
    A1(12:211,12:211)=W;
    N=zeros(222);
    N(:)=100;
    N(12:211,12:211)=M;
    for x=12:211
        for y=12:211
            if A1(x,y)<100
            for i=1:11
                C=zeros(23);
                C(12-i+1:12+i-1:12-i+1:12+i-1)=1;
                B=C(12-i:12+i,12-i:12+i);
                if(A1(x-i:x+i,y-i:y+i))==100|B
                    [p,q]=find(A1(x-i:x+i:y-i:y+i)<100);
                    p=p-1+x-i;
                    q=q-1+y-i;
                    m=fix(sum(p)/length(p));
                    n=fix(sum(q)/length(q));   %m，n為重心坐標(biāo)
                    G=N(x-i:x+i,y-i:y+i);
                    D(m,n)=100*(2*i+1)*(2*i+1)-sum(G(:));%重心上的粒子數(shù)
                    E(m,n)=A1(x,y)-length(p);%重心上的能量
                    if E(m,n)<0
                        E(m,n)=0;
                    end
                    break
                end
            end
        end
    end
end