Vortex pair.py

From Werner KRAUTH

This is the program Vortex_pair.py, that you need for Homework 10 of the 2019 ICFP lectures. Don't hesitate to copy, paste and modify this program that was originally written by Ze Lei.

#
# Based on an original program written by Ze LEI (LPS ENS)
#
import matplotlib.pyplot as plt
import numpy
import math

def plot_spin(data, L):
    data = data.reshape(L*L)
    x = numpy.arange(0, L, 1)
    y = numpy.arange(0, L, 1)
    X, Y = numpy.meshgrid(x, y)
    U = numpy.cos(data)
    V = numpy.sin(data)
    data = data.reshape(L, L)
    U = U.reshape(L, L)
    V = V.reshape(L, L)

    width = 0.5
    plt.figure()
    plt.quiver(X, Y, U, V, data, pivot='mid', headlength=10, linewidths=width)
    plt.title('Vortex_antivortex_L=' + str(L) + '_offset=' +  str(offset))
    plt.savefig('Vortex_antivortex_L=' + str(L) + '_offset=' +  str(offset)  + '.pdf')
    plt.close()

def xy_to_phi(xx, yy):
    phi = math.atan(yy / xx)
    if xx < 0.0: phi += math.pi
    return phi

L = 64 # use even values of L

XX = numpy.zeros((L, L)); YY = numpy.zeros((L, L))
phi = numpy.zeros((L, L))

#vortex position
offset = 6 # attention, "offset" is not the cartesian distance between + and -
core_v = numpy.array([L / 2 - offset, L / 2 - offset]) - 0.5
core_a = numpy.array([L / 2 + offset, L / 2 + offset]) - 0.5

#initial configuration with vortex
for xx in range(L):
    for yy in range(L):
        #vortex
        Xv =  -(yy - core_v[0])
        Yv =  -(xx - core_v[1])
        Nv = Xv*Xv + Yv*Yv
        #antivortex
        Xa = xx - core_a[0]
        Ya = yy - core_a[1]
        Na = Xa*Xa + Ya*Ya
        Xx = Xv / Nv + Xa / Na
        Yy = Yv / Nv + Ya / Na
        Norm = math.sqrt(Xx*Xx + Yy*Yy)
        XX[xx, yy] = Xx/Norm; YY[xx, yy] = Yy/Norm

#Neighbor relations:
neighbours = [[[] for xx in range(L)] for yy in range(L)]
for xx in range(L):
    for yy in range(L):
        if xx > 0: neighbours[xx][yy].append((xx-1, yy))
        if xx < L - 1: neighbours[xx][yy].append((xx+1, yy))
        if yy > 0: neighbours[xx][yy].append((xx, yy-1))
        if yy < L - 1: neighbours[xx][yy].append((xx, yy+1))

#Minimization:
for i in range(0, 1000):
    for xx in range(L):
        for yy in range(L):
            Xx = sum([XX[site] for site in neighbours[xx][yy] ])
            Yy = sum([YY[site] for site in neighbours[xx][yy] ])
            Norm = math.sqrt(Xx * Xx + Yy * Yy)
            XX[xx, yy] = Xx / Norm
            YY[xx, yy] = Yy / Norm
    if i % 10 == 0:
        energy = 2.0 * L * (L - 1) - numpy.sum(XX[:, 0:L - 1] * XX[:, 1:L] + YY[:,
              0:L-1] * YY[:, 1:L]) - numpy.sum(XX[0:L-1, :] * XX[1:L, :] + YY[0:L-1, :] * YY[1:L, :])
        print i, energy
for xx in range(L):
    for yy in range(L):
        phi[xx, yy] = xy_to_phi(XX[xx, yy], YY[xx, yy])
        phi[xx, yy] = (phi[xx, yy] + 2 * math.pi) % (2 * math.pi)
plot_spin(phi, L)