function 隣接格子点(ix,iy,Lx,Ly,周期境界)
    const lmax = 4 #最近接格子点の数
    ls = Array{Array{Int64,1}}(lmax)
    ds = ((+1,0),(-1,0),(0,+1),(0,-1))
    count = 0
    for d in ds        
        jx = ix + d[1]
        jy = iy + d[2]
        
        if 周期境界
            jx += ifelse(jx>Lx,-Lx,0)
            jx += ifelse(jx<1,Lx,0)    
            jy += ifelse(jy>Ly,-Ly,0)
            jy += ifelse(jy<1,Ly,0)                                  
        end     
        
        if 1 <= jx <= Lx && 1 <= jy <= Ly
            count += 1
            ls[count] = [jx,jy]
        end        
    end    
    return ls[1:count]

end

function 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
    #ls = 隣接格子点(ix,iy,Lx,Ly,周期境界)
    S = 0
    if ix == 1
        left = スピン配置[Lx,iy]
    else
        left = スピン配置[ix-1,iy]
    end
    if ix == Lx
        right = スピン配置[1,iy]
    else
        right = スピン配置[ix+1,iy]
    end
    if iy == 1
        down = スピン配置[ix,Ly]
    else
        down = スピン配置[ix,iy-1]
    end
    if iy == Ly
        up = スピン配置[ix,1]
    else
        up = スピン配置[ix,iy+1]
    end

    S = (left+right+down+up)

    return S
end

function 全エネルギー(スピン配置,周期境界,J,h)
    Lx = size(スピン配置,1)
    Ly = size(スピン配置,2)        
    E = 0.0
    for ix in 1:Lx
        for iy in 1:Ly
            σi = スピン配置[ix,iy]
            S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)            
            E += -(J/2)*σi*S-h*σi
        end
    end
    return E
end

function メトロポリス法!(σi,ΔE,T,E,mz)
    σ_new = ifelse(rand() <= exp(-ΔE/T),-σi,σi)
    E,mz,accept = ifelse(σ_new == -σi,(E+ΔE,mz-2*σi,1),(E,mz,0))
    return σ_new,E,accept,mz
end

function 熱浴法!(σi,ΔE,T,E,mz)
    α = σi*ΔE
    ratio = (1 + tanh(α/2T))/2 #アップスピンの受け入れ確率    
    σ_new = ifelse(rand() <= ratio,+1,-1)
    E,mz,accept = ifelse(σ_new == -σi,(ΔE,-2*σi,1),(E,mz,0))
    return σ_new,E,accept,mz    
end

function 配置アップデート!(スピン配置,周期境界,flipmethod,T,J,h,E,mz)
    Lx = size(スピン配置,1)
    Ly = size(スピン配置,2)    
    ix = rand(1:Lx)
    iy = rand(1:Ly)  
    σi = スピン配置[ix,iy]
    S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
    ΔE =  2J*σi*S+2h*σi
    スピン配置[ix,iy],E,accept,mz = flipmethod(σi,ΔE,T,E,mz)
    #println(accept)
    return スピン配置,E,mz,accept
end

using ProgressMeter
function モンテカルロ法!(スピン配置,周期境界,flipmethod,T,J,h,最大ステップ,熱化ステップ)
    prog = Progress(最大ステップ,1)
    E = 全エネルギー(スピン配置,周期境界,J,h)
    mz = sum(スピン配置)
    accept = 0
    受けいれ総数 = 0
    totalcount = 最大ステップ-熱化ステップ
    Mz = Array{Float64}(totalcount)
    Energy = Array{Float64}(totalcount)
    比熱 = Array{Float64}(totalcount)
    受けいれ確率 = Array{Float64}(totalcount)
    Mz2 = Array{Float64}(totalcount) #磁化の二乗
    mzヒストグラム = zeros(Int64,length(スピン配置)*2+1)
    mzsum = 0.0
    mz2sum = 0.0
    energy = 0.0
    energy2 = 0.0
    count = 0
    for i=1:最大ステップ
        for j=1:length(スピン配置)
            スピン配置,E,mz,accept = 配置アップデート!(スピン配置,周期境界,flipmethod,T,J,h,E,mz) 
        end
        if rand() < 0.01            
            スピン配置 = -スピン配置
            mz = -mz
            accept = 1
        end
        if i > 熱化ステップ 
            mzヒストグラム[mz+length(スピン配置)+1] += 1
            count += 1 
            mzsum += mz/length(スピン配置) 
            mz2sum += (mz/length(スピン配置))^2             
            energy += E
            energy2 += E^2
            受けいれ総数 += accept 
            Mz[count] = mzsum/count
            Mz2[count] = mz2sum/count            
            Energy[count] = energy/count
            比熱[count] = (energy2/count-Energy[count]^2)/T^2
            受けいれ確率[count] = 受けいれ総数/count
        end 
        next!(prog)
    end
    return Mz[1:count],Mz2[1:count],Energy[1:count],比熱[1:count],受けいれ確率[1:count],mzヒストグラム

end

using Plots





Tc = 2/log(1+sqrt(2))

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 1.0   
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

function 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
    S = スピン配置[(ix-2+Lx)%Lx+1,iy]+スピン配置[ix%Lx+1,iy]+スピン配置[ix,(iy-2+Ly)%Ly+1]+スピン配置[ix,iy%Ly+1]
    return S
end

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 1.0   
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

function 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
    S = 0
    jx = ifelse(ix==Lx,1,ix+1)
    jy = iy
    S += スピン配置[jx,jy]
    jx = ifelse(ix==1,Lx,ix-1)
    jy = iy 
    S += スピン配置[jx,jy]
    jy = ifelse(iy==Ly,1,iy+1)
    jx = ix
    S += スピン配置[jx,jy]
    jy = ifelse(iy==1,Ly,iy-1)
    jx = ix 
    S += スピン配置[jx,jy]    
    
    
#    S = スピン配置[(ix-2+Lx)%Lx+1,iy]+スピン配置[ix%Lx+1,iy]+スピン配置[ix,(iy-2+Ly)%Ly+1]+スピン配置[ix,iy%Ly+1]
    return S
end

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 1.0   
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,J,h,E,mz)
    accept = 0
    acc = 0
    for j=1:Lx*Ly
        スピン配置,E,mz,acc = 配置アップデート!(スピン配置,周期境界,flipmethod,T,J,h,E,mz)
        accept += acc
    end
    return スピン配置,E,mz,accept
end

using ProgressMeter
function モンテカルロ法!(スピン配置,周期境界,flipmethod,T,J,h,最大ステップ,熱化ステップ)
    prog = Progress(最大ステップ,1)
    E = 全エネルギー(スピン配置,周期境界,J,h)
    mz = sum(スピン配置)
    accept = 0
    受けいれ総数 = 0
    totalcount = 最大ステップ-熱化ステップ
    Mz = Array{Float64}(totalcount)
    Energy = Array{Float64}(totalcount)
    比熱 = Array{Float64}(totalcount)
    受けいれ確率 = Array{Float64}(totalcount)
    Mz2 = Array{Float64}(totalcount) #磁化の二乗
    mzヒストグラム = zeros(Int64,length(スピン配置)*2+1)
    mzsum = 0.0
    mz2sum = 0.0
    energy = 0.0
    energy2 = 0.0
    count = 0
    Lx = size(スピン配置,1)
    Ly = size(スピン配置,2)
    for i=1:最大ステップ
        スピン配置,E,mz,accept = スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,J,h,E,mz)
        if rand() < 0.01            
            スピン配置 = -スピン配置
            mz = -mz
            accept = 1
        end
        if i > 熱化ステップ 
            mzヒストグラム[mz+length(スピン配置)+1] += 1
            count += 1 
            mzsum += mz/length(スピン配置) 
            mz2sum += (mz/length(スピン配置))^2             
            energy += E
            energy2 += E^2
            受けいれ総数 += accept 
            Mz[count] = mzsum/count
            Mz2[count] = mz2sum/count            
            Energy[count] = energy/count
            比熱[count] = (energy2/count-Energy[count]^2)/T^2
            受けいれ確率[count] = 受けいれ総数/count
        end 
        next!(prog)
    end
    return Mz[1:count],Mz2[1:count],Energy[1:count],比熱[1:count],受けいれ確率[1:count],mzヒストグラム

end

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 1.0   
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,J,h,E,mz)
    accept = 0
    acc = 0
    #for j=1:Lx*Ly
    for jx in 1:Lx
        for jy in 1:Ly
            ix = rand(1:Lx)
            iy = rand(1:Ly)          
            σi = スピン配置[ix,iy]    
            S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
            ΔE =  2J*σi*S+2h*σi
            σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz)            
            スピン配置[ix,iy] = σ_new
            accept += acc
        end
    end
    return スピン配置,E,mz,accept
end


Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 1.0   
J = 1.0
h = 0.0
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,J,h,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

function boltzmann(T)
    w = zeros(Float64,8*2+1)    
    for ΔE in -8:8
        w[ΔE+9] = exp(-ΔE/T)
    end
    return w
end

function メトロポリス法!(σi,ΔE,T,E,mz,w)
    σ_new = ifelse(rand() <= w[ΔE+9],-σi,σi)
    E,mz,accept = ifelse(σ_new == -σi,(E+ΔE,mz-2*σi,1),(E,mz,0))
    return σ_new,E,accept,mz
end

function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w) #磁場hとカップリングJを除去。wを追加
    accept = 0
    acc = 0
    #for j=1:Lx*Ly
    for jx in 1:Lx
        for jy in 1:Ly
            ix = rand(1:Lx)
            iy = rand(1:Ly)          
            σi = スピン配置[ix,iy]    
            S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
            ΔE =  2σi*S
            σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz,w)#wを追加            
            スピン配置[ix,iy] = σ_new
            accept += acc
        end
    end
    return スピン配置,E,mz,accept
end

function モンテカルロ法!(スピン配置,周期境界,flipmethod,T,最大ステップ,熱化ステップ)
    prog = Progress(最大ステップ,1)
    w = boltzmann(T) #追加した部分
    E = 全エネルギー(スピン配置,周期境界,J,h)
    mz = sum(スピン配置)
    accept = 0
    受けいれ総数 = 0
    totalcount = 最大ステップ-熱化ステップ
    Mz = Array{Float64}(totalcount)
    Energy = Array{Float64}(totalcount)
    比熱 = Array{Float64}(totalcount)
    受けいれ確率 = Array{Float64}(totalcount)
    Mz2 = Array{Float64}(totalcount) #磁化の二乗
    mzヒストグラム = zeros(Int64,length(スピン配置)*2+1)
    mzsum = 0.0
    mz2sum = 0.0
    energy = 0.0
    energy2 = 0.0
    count = 0
    Lx = size(スピン配置,1)
    Ly = size(スピン配置,2)
    for i=1:最大ステップ
        スピン配置,E,mz,accept = スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w)
        if rand() < 0.01            
            スピン配置 = -スピン配置
            mz = -mz
            accept = 1
        end
        if i > 熱化ステップ 
            mzヒストグラム[mz+length(スピン配置)+1] += 1
            count += 1 
            mzsum += mz/length(スピン配置) 
            mz2sum += (mz/length(スピン配置))^2             
            energy += E
            energy2 += E^2
            受けいれ総数 += accept 
            Mz[count] = mzsum/count
            Mz2[count] = mz2sum/count            
            Energy[count] = energy/count
            比熱[count] = (energy2/count-Energy[count]^2)/T^2
            受けいれ確率[count] = 受けいれ総数/count
        end 
        next!(prog)
    end
    return Mz[1:count],Mz2[1:count],Energy[1:count],比熱[1:count],受けいれ確率[1:count],mzヒストグラム

end

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 1.0   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = 3.0   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000
熱化ステップ = 200

T = Tc   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
熱化ステップ = 200

T = Tc   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

function メトロポリス法!(σi,ΔE,T,E,mz,w,r)
    σ_new = ifelse(r <= w[ΔE+9],-σi,σi)
    E,mz,accept = ifelse(σ_new == -σi,(E+ΔE,mz-2*σi,1),(E,mz,0))
    return σ_new,E,accept,mz
end


function スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w) #磁場hとカップリングJを除去。wを追加
    r = rand(Lx,Ly)    
    accept = 0
    acc = 0
    #for j=1:Lx*Ly
    for jx in 1:Lx
        for jy in 1:Ly
            ix = jx#rand(1:Lx)
            iy = jy#rand(1:Ly)          
            σi = スピン配置[ix,iy]    
            S = 周辺スピン和(ix,iy,Lx,Ly,スピン配置,周期境界)
            ΔE =  2σi*S
            #σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz,w)#wを追加    
            σ_new,E,acc,mz = flipmethod(σi,ΔE,T,E,mz,w,r[ix,iy])#wを追加             
            スピン配置[ix,iy] = σ_new
            accept += acc
        end
    end
    return スピン配置,E,mz,accept
end



Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
熱化ステップ = 200

T = Tc   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 1000000
熱化ステップ = 200

T = Tc   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time = Int64[]
for i in 1:length(mzヒストグラム)
    push!(time,i-Lx*Ly-1)
end 
plot(time./(Lx*Ly),mzヒストグラム)

Lx = 50
Ly = 50
最大ステップ = 1000000
熱化ステップ = 200

srand(123)
スピン配置 = rand([-1,1],Lx,Ly)

T = Tc   
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム1 = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
time1 = Int64[]
mzヒストグラム1 = length(mzヒストグラム1)*mzヒストグラム1/sum(mzヒストグラム1)
for i in 1:length(mzヒストグラム1)
    push!(time1,i-Lx*Ly-1)
end 
time1 = time1./(Lx*Ly)
plot(time1,mzヒストグラム1)

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
 
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム2 = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
mzヒストグラム2 = length(mzヒストグラム2)*mzヒストグラム2/sum(mzヒストグラム2)
time2 = Int64[]
for i in 1:length(mzヒストグラム2)
    push!(time2,i-Lx*Ly-1)
end 
time2 = time2./(Lx*Ly)
plot!(time2,mzヒストグラム2)

Lx = 150
Ly = 150
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)

  
@time Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム3 = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
mzヒストグラム3 = length(mzヒストグラム3)*mzヒストグラム3/sum(mzヒストグラム3)
time3 = Int64[]
for i in 1:length(mzヒストグラム3)
    push!(time3,i-Lx*Ly-1)
end 
time3 = time3./(Lx*Ly)
plot!(time3,mzヒストグラム3)

Lx = 50
Ly = Lx
plot((time1)*(Lx^(1/8)),mzヒストグラム1,title="1000000 steps",label=string(Lx)*"x"*string(Ly))
Lx = 100
Ly = Lx
plot!((time2)*(Lx^(1/8)),mzヒストグラム2,title="1000000 steps",label=string(Lx)*"x"*string(Ly))
Lx = 150
Ly = Lx
plot!((time3)*(Lx^(1/8)),mzヒストグラム3,title="1000000 steps",label=string(Lx)*"x"*string(Ly))

function モンテカルロ法可視化!(スピン配置,周期境界,flipmethod,T,最大ステップ,filename)
    E = 全エネルギー(スピン配置,周期境界,J,h)
    w = boltzmann(T) #追加した部分
    mz = sum(スピン配置)
    accept = 0
    prog = Progress(最大ステップ,1)
    
    maps = @animate for i=1:最大ステップ
        スピン配置,E,accept,mz = スイープ!(スピン配置,Lx,Ly,周期境界,flipmethod,T,E,mz,w)
        heatmap(1:Lx, 1:Ly, スピン配置,aspect_ratio=:equal)
        next!(prog)
        end every 100
    gif(maps, "./"*filename, fps = 15)   

end

Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000

T = 1.0   
J = 1.0
filename = "ising_fps15T1.gif"
@time  モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)



Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000

T = Tc   
J = 1.0
filename = "ising_fps15Tc.gif"
@time  モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)



Lx = 100
Ly = 100
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 10000

T = 3.0   
J = 1.0
filename = "ising_fps15T3.gif"
@time  モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)



Lx = 200
Ly = 200
周期境界 = true
srand(123)
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000

T = 1.0   
J = 1.0
filename = "ising_fps15T1_200.gif"
@time  モンテカルロ法可視化!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,filename)



srand(123)
Lx = 50
Ly = 50
スピン配置 = rand([-1,1],Lx,Ly)
最大ステップ = 100000
熱化ステップ = 200
Mztdep = Float64[]
比熱tdep = Float64[]
tdep = Float64[]
@time for i in 1:20
    T = (21-i)*0.2
    Mz,Mz2,Energy,比熱,受けいれ確率,mzヒストグラム = モンテカルロ法!(スピン配置,周期境界,メトロポリス法!,T,最大ステップ,熱化ステップ)
    push!(Mztdep,sqrt(Mz2[end]))
    push!(比熱tdep,比熱[end])    
    push!(tdep,T)
end
plot(tdep,Mztdep)

plot(tdep,比熱tdep)

filename = "./dc.dat"
fd = open( filename, "r" )
温度f = zeros(Float64,20)
比熱f = zeros(Float64,20)
磁化f = zeros(Float64,20)
cnt = 0
while !eof(fd) 
    cnt += 1
    line = readline(fd)
    if cnt >3 
        u = split(line)

        温度f[cnt-3] = parse(Float64,u[1])
        比熱f[cnt-3] = parse(Float64,u[3])
        磁化f[cnt-3] = parse(Float64,u[7])
    end
end

plot(温度f,磁化f,marker=:circle,label="Fortran")
plot!(tdep,Mztdep,marker=:circle,label="Julia")