Optimizing my 2D Ising model code in Julia I'm just starting to learn Julia, I work primarily in physics and am used to writing most of my code in Fortran90 and occasionally Python for Tensorflow (also Mathematica but that's less relevant). Julia has been recommended to me and I started checking it out; I like it a lot in theory as a middle ground between the speed of Fortran and the syntax of Python. To test it out I wrote a simple 2D Ising model code implementing a basic single-spin-flip Metropolis Monte Carlo algorithm. However, this code runs very slowly compared to an equivalent code in Fortran. Am I doing something wrong which is significantly affecting the performance of the code? I know almost nothing beyond what I've done here. I am using the Juno IDE in Atom on Windows 10. As an aside, I would also like to know how I can make multiple plots in the Atom plot tab, but that's secondary.

using Printf
using Plots

L       = 20             # linear size of lattice
n_sweep = 20             # number of sweeps between sampling
n_therm = 1000           # number of sweeps to thermalize
n_data  = 100            # number of data samples per temperature
temps   = 4.0:-0.3:0.1   # temperatures to sample
e1 = Array(1:n_data)     # array to hold energy measurements (fixed T)
m1 = Array(1:n_data)     # array to hold magnetization measurements (fixed T)
et = []                  # array to append average energy at each T
mt = []                  # "                      magnetizations
s  = ones(Int32,L,L)     # lattice of Ising spins (+/-1)

function measure(i)      # measure i'th sample of energy and magnetization
    en = 0
    m = 0
    for x = 1:L
        for y = 1:L
            u = 1+mod(y,L) # up
            r = 1+mod(x,L) # right
            en -= s[x,y]*(s[x,u]+s[r,y]) # energy
            m  += s[x,y]                 # magnetization
        end
    end
    energy[i] = en
    magnetization[i] = abs(m)
end

function flip(x,y,T) # apply metropolis spin flip algorithm to site (x,y) w/ temp T
    u = 1+mod(y,L)   # up
    d = 1+mod(y-2,L) # down
    r = 1+mod(x,L)   # right
    l = 1+mod(x-2,L) # left
    de = 2*s[x,y]*(s[x,u]+s[x,d]+s[l,y]+s[r,y])
    if (de &lt; 0)
        s[x,y] = -s[x,y]
    else
        p = rand()
        if (p &lt; exp(-de/T))
            s[x,y] = -s[x,y]
        end
    end
end

function sweep(n,T) # apply flip() to every site on the lattice
    for i = 1:n
        for x = 1:L
            for y = 1:L
                flip(x,y,T)
            end
        end
    end
end

for T in temps              # loop over temperatures
    sweep(n_therm, T)       # thermalize the lattice
    energy        = e1      # reset energy measurement array
    magnetization = m1      # same
    for i = 1:n_data        # take n_data measurements w/ n_sweep 
        sweep(n_sweep, T)   
        measure(i)
    end
    en_ave = sum(energy)/n_data           # compute average
    ma_ave = sum(magnetization)/n_data
    push!(et,en_ave/(L*L))                # add to the list
    push!(mt,ma_ave/(L*L))
    @printf("%8.3f  %8.3f \n", en_ave/(L*L), ma_ave/(L*L))
end

plot(temps,mt) # plot magnetization vs. temperature
#plot(temps,et)
The biggest problem with this code is the amount of things in global scope. Rewriting it to make some things const, and passing in the rest from a main function brings the time down to .7 seconds (not including the using plots which takes 3.5 seconds, compared to ~10 seconds before. Here is the updated code, I hope it helps.

using Printf
using Plots

const L       = 20             # linear size of lattice
const n_sweep = 20             # number of sweeps between sampling
const n_therm = 1000           # number of sweeps to thermalize
const n_data  = 100            # number of data samples per temperature
const temps   = 4.0:-0.3:0.1   # temperatures to sample

function measure(i, energy, magnetization, s)      # measure i'th sample of energy and magnetization
    en = 0
    m = 0
    for x = 1:L
        for y = 1:L
            u = 1+mod(y,L) # up
            r = 1+mod(x,L) # right
            en -= s[x,y]*(s[x,u]+s[r,y]) # energy
            m  += s[x,y]                 # magnetization
        end
    end
    energy[i] = en
    magnetization[i] = abs(m)
end

function flip(x, y, T, s) # apply metropolis spin flip algorithm to site (x,y) w/ temp T
    u = 1+mod(y,L)   # up
    d = 1+mod(y-2,L) # down
    r = 1+mod(x,L)   # right
    l = 1+mod(x-2,L) # left
    de = 2*s[x,y]*(s[x,u]+s[x,d]+s[l,y]+s[r,y])
    if (de &lt; 0)
        s[x,y] = -s[x,y]
    else
        p = rand()
        if (p &lt; exp(-de/T))
            s[x,y] = -s[x,y]
        end
    end
end

function sweep(n, T, s) # apply flip() to every site on the lattice
    for i = 1:n
        for x = 1:L
            for y = 1:L
                flip(x,y,T, s)
            end
        end
    end
end

function main()
    e1 = Array(1:n_data)     # array to hold energy measurements (fixed T)
    m1 = Array(1:n_data)     # array to hold magnetization measurements (fixed T)
    et = []                  # array to append average energy at each T
    mt = []                  # "                      magnetizations
    s  = ones(Int32,L,L)     # lattice of Ising spins (+/-1)
    for T in temps              # loop over temperatures
        sweep(n_therm, T, s)    # thermalize the lattice
        energy        = e1      # reset energy measurement array
        magnetization = m1      # same
        for i = 1:n_data        # take n_data measurements w/ n_sweep 
            sweep(n_sweep, T, s)   
            measure(i, energy, magnetization, s)
        end
        en_ave = sum(energy)/n_data           # compute average
        ma_ave = sum(magnetization)/n_data
        push!(et,en_ave/(L*L))                # add to the list
        push!(mt,ma_ave/(L*L))
        @printf("%8.3f  %8.3f \n", en_ave/(L*L), ma_ave/(L*L))
    end
    plot(temps,mt) # plot magnetization vs. temperature
    #plot(temps,et)
end

main()