# Vortex pair.py

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)
```