% Ising模型估計U和C
function [U,C]=Ising(M,p,N,N1,N2,number,preselect,image)
%函數(shù)格式:[U,C] =  Ising(n,p,N,N1,N2,times,preselect,image)
%輸入:    n - 二維Ising模型的邊長,滿足周期邊界條件
%        p - beta的取值個數(shù)
%        N,N1,N2,times - 組策略,具體地,N為馬氏鏈長度,N1為求beta中舍去的開頭的
%        狀態(tài)數(shù),N2為取樣本(狀態(tài))的間隔,times為計算次數(shù)
%        preselect - 預(yù)選方式,1=等可能預(yù)選,2=單點變化預(yù)選
%        image - 選擇輸出輸出的圖像,1=輸出圖像,其它,不輸出
%輸出:   U - 內(nèi)能U的結(jié)果序列
%        C - 比熱容C的結(jié)果序列

n=M-1;
U([1:p+1])=0;
C([1:p+1])=0;

for k=1:number
    for i=0:p
        beta=i/p;
        sigma=initiate(n);
        [H,H2]=FuncHSpace(sigma,N,N1,N2,beta,preselect);
        U(i+1)=U(i+1)+H/(n+1)/(n+1);
        C(i+1)=C(i+1)+beta^2/(n+1)^2*(H2-H*H);
    end
end

U=U/number;
C=C/number;

if image==1
    figure;
    subplot(1,2,1);
    plot([0:1/p:1],U,'r','linewidth',2);
    xlabel('\beta');
    ylabel('U');
    %figure;
    subplot(1,2,2);
    plot([0:1/p:1],C,'b','linewidth',2);
    xlabel('\beta');
    
    ylabel('C');
end

end