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 < 0) s[x,y] = -s[x,y] else p = rand() if (p < 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)