versioninfo()

function スイープ!(配置,Lx,Ly,w,mz)
    for ix in 1:Lx
        for iy in 1:Ly  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
        end
    end
    return 配置,mz
end

function モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
    w = exp.(-[i for i in -8:8]/T) #ボルツマン重み e^[-ΔE/T]
    mz = sum(配置)
    Mz = 0
    for i=1:熱化ステップ
        配置,mz=スイープ!(配置,Lx,Ly,w,mz)        
    end    
    for i=1:最大ステップ
        配置,mz=スイープ!(配置,Lx,Ly,w,mz)
        Mz += abs(mz)
    end
    return Mz/(Lx*Ly*最大ステップ)
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

function スイープ!(配置,Lx,Ly,w,mz)
    r = rand(Lx,Ly) #乱数確保
    for ix in 1:Lx
        for iy in 1:Ly  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
            配置[ix,iy] = ifelse(r[ix,iy] <= w[ΔE+9],-σi,σi)
            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
        end
    end
    return 配置,mz
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

function スイープ!(配置,Lx,Ly,w,mz)
    for ix in 1:Lx
        for iy in 1:Ly  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ix==Lx?1:ix+1,iy]+配置[ix==1?Lx:ix-1,iy]+配置[ix,iy==Ly?1:iy+1]+配置[ix,iy==1?Ly:iy-1])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
        end
    end
    return 配置,mz
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

function スイープ!(配置,Lx,Ly,w,mz)
    for ix in 1:Lx
        for iy in 1:Ly  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
        end
    end
    return 配置,mz
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

ix = 1
println(mod1(ix-1,Lx))
ix = Lx
println(mod1(ix+1,Lx))

function スイープ!(配置,Lx,Ly,w,mz)
    for ix in 1:Lx
        for iy in 1:Ly  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[mod1(ix+1,Lx),iy]+配置[mod1(ix-1,Lx),iy]+配置[ix,mod1(iy+1,Ly)]+配置[ix,mod1(iy-1,Ly)])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
        end
    end
    return 配置,mz
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

function スイープ!(配置,Lx,Ly,w,mz)
    ix = 1
    iy = 1
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,Ly])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  
    
    for iy in 2:Ly-1
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,iy-1])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)        
    end
    iy = Ly
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,1]+配置[ix,iy-1])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)      
        
    for ix in 2:Lx-1
        iy = 1
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  
        
        for iy in 2:Ly-1  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ix+1,iy]+配置[iy-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
        end
        
        iy = Ly
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)             
    end

    ix = Lx
    iy = 1
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  
    
    for iy in 2:Ly-1
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)        
    end
    iy = Ly
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  
    
    
    return 配置,mz
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

function スイープ!(配置,Lx,Ly,w,mz)
    ix = 1
    iy = 1
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,Ly])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += 配置[ix,iy]-σi
    
    for iy in 2:Ly-1
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,iy-1])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += 配置[ix,iy]-σi       
    end
    iy = Ly
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,1]+配置[ix,iy-1])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += 配置[ix,iy]-σi     
        
    for ix in 2:Lx-1
        iy = 1
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += 配置[ix,iy]-σi
        
        for iy in 2:Ly-1  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ix+1,iy]+配置[iy-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += 配置[ix,iy]-σi
        end
        
        iy = Ly
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += 配置[ix,iy]-σi          
    end

    ix = Lx
    iy = 1
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += 配置[ix,iy]-σi 
    
    for iy in 2:Ly-1
        σi = 配置[ix,iy]
        ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
        mz += 配置[ix,iy]-σi      
    end
    iy = Ly
    σi = 配置[ix,iy]
    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
    mz += 配置[ix,iy]-σi
    
    
    return 配置,mz
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

function スイープ!(配置,Lx,Ly,w,mz)
    for ix in 1:Lx
        for iy in 1:Ly  
            σi = 配置[ix,iy]
            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
            mz += 配置[ix,iy]-σi
        end
    end
    return 配置,mz
end


function モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
    w = exp.(-[i for i in -8:8]/T) #ボルツマン重み e^[-ΔE/T]
    mz = sum(配置)
    Mz = 0
    for i=1:熱化ステップ
        配置,mz=スイープ!(配置,Lx,Ly,w,mz)        
    end    
    for i=1:最大ステップ
        配置,mz=スイープ!(配置,Lx,Ly,w,mz)
        Mz += abs(mz)
    end
    return Mz/(Lx*Ly*最大ステップ)
end

Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)

using Plots
pyplot()

Lx = 100
Ly = 100
最大ステップ = 100000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
nt = 20
time = [(nt+1-i)*0.2 for i in 1:nt]
Mztdep = zeros(Float64,nt)
i = 0
@time for T in time
    i += 1
    Mztdep[i]  = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
end

plot(time,Mztdep)