---
jupytext:
  formats: md:myst,ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
import mmf_setup;mmf_setup.nbinit()
```

Here we derive some numerical techniques for imprinting a vortex dipole in a periodic box after [Billam:2014].

[Billam:2014]: https://doi.org/10.1103/physrevlett.112.145301 'T. P. Billam, M. T. Reeves, B. P. Anderson, and A. S. Bradley, "Onsager-Kraichnan Condensation in Decaying Two-Dimensional Quantum Turbulence", Phys. Rev. Lett. 112(14),  (2014) '

+++

# The Problem

+++

In 2D infinite matter, a vortex at position $z_0 = x_0 + \I y_0$ can be described, up to an overall phase, by

$$
  \psi(z) = f(\abs{z-z_0}) (z-z_0) \propto e^{\I\phi(z)}
$$

where $f(r) \sim r^2$ for small $r \ll \xi$ describes the core structure.  For our purposes, we will neglect this dependence and focus on the form of the phase $\phi(z)$:

$$
  \phi(z) = \tan^{-1}\frac{y-y_0}{x-x_0}.
$$

This structure is generally incompatible with periodic boundary conditions due to the global rotation (but see...).

If we consider istead a vortex dipole with an additional anti-vortex at $z_1 = x_1 + \I y_1$, then we can produce a periodic lattice by considering all images:

$$
  \psi_{D}(z) \propto  (z-z_0)(z-z_1)^{*}, \qquad
  \psi(z) \propto \prod_{n_x, n_y} \psi_{D}\Bigl(z-L(n_x,n_y)\Bigr), \qquad
  L(n_x, n_y) = n_xL_x + \I n_yL_y.
$$

```{code-cell}
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from tqdm import tqdm
from itertools import product
Lx = 2*np.pi
Ly = 2*np.pi
d = 2.0
z0 = Lx/4 + 0.25j*Ly
z1 = Lx/4+d + (0.25*Ly+d)*1j

def psi_D(z):
    psi = (z-z0)*(z-z1).conjugate()
    return psi/abs(psi)

def get_psi(z, Nx=12, Ny=12):
    psi = np.prod([
        psi_D(z-(nx*Lx+1j*ny*Ly)) 
        for nx, ny in (product(
            range(-Nx//2, Nx//2+1),
            range(-Ny//2, Ny//2+1)))], axis=0)
    return psi/abs(psi)
```

```{code-cell}
N = 100
abs(get_psi(np.pi + 2j*np.pi, Nx=N) - get_psi(np.pi + 0j*np.pi, Nx=N))
```

```{code-cell}
x = np.linspace(0, Lx, 100)
y = np.linspace(0, Ly, 102)
z = x[:, None] + 1j*y[None, :]
psi = get_psi(z, Nx=N, Ny=N)
plt.pcolormesh(x, y, np.angle(psi).T)
```

```{code-cell}
Ns = 2**np.arange(2, 8)
psis = np.array([get_psi(z, Nx=N, Ny=N) for N in Ns])
plt.loglog(1/Ns, abs(psis - psis[-1]).max(axis=(-1)).max(axis=-1), '+')
```

$$
  \psi = \psi_0 + a\frac{1}{N^2}
$$

```{code-cell}
shape = (psis.shape[0], )
N = 5
psi0 = np.polyfit(1/Ns[:N]**2, psis[:N].reshape((N, np.prod(psis.shape[1:]))), 4)[-1].reshape(psis.shape[1:])
plt.plot(np.log(Ns), np.log(abs(psis-psi0).max(axis=(-1)).max(axis=-1)), '+')
plt.plot(np.log(Ns), -2.3-2*np.log(Ns))
```

```{code-cell}
Ns = 2**np.arange(2, 9)
psis = np.array([get_psi(z, Nx=N, Ny=N) for N in Ns])
plt.loglog(1/Ns, abs(psis - psi0).max(axis=(-1)).max(axis=-1), '+')
```

```{code-cell}
def A(x):
    return np.arctan(x)

_x = 0.1
_a = 0.2
sum(A(_x/((2*k-1)*np.pi/2 - _a)) - A(_x/((2*k-1)*np.pi/2 + _a)) for k in range(1, 100)), A(np.tanh(_x)*np.tan(_a))
```

```{code-cell}
x, y = z.real, z.imag
x0, y0, x1, y1 = z0.real, z0.imag, z1.real, z1.imag
X0, Y0, X1, Y1 = x-x0, y-y0, x-x1, y-y1


Eq5 = sum(A(Y1/((2*m - 1)*np.pi + (X1 - np.pi))) 
        - A(Y1/((2*m - 1)*np.pi - (X1 - np.pi)))
          for m in range(1, 1000))
Eq6 = -A(np.tanh(Y1/2) * np.tan((X1 - np.pi)/2))
plt.imshow(Eq6-Eq5)
plt.colorbar()
```

```{code-cell}
# Formula from [Billam:2014]
def get_theta(z, z0, z1, N=5):
    # X^+ = X0, X^- = X1
    x, y = z.real, z.imag
    x0, y0, x1, y1 = z0.real, z0.imag, z1.real, z1.imag    
    X0, Y0, X1, Y1 = x-x0, y-y0, x-x1, y-y1
    res = (
        sum(np.arctan(np.tanh((Y1 + 2*np.pi*n)/2) * np.tan((X1 - np.pi)/2))
            - np.arctan(np.tanh((Y0 + 2*np.pi*n)/2) * np.tan((X0 - np.pi)/2))
            + np.pi * (1.0*(X0>0) - (X1>0))
            for n in range(-N, N+1))
        - (x0 - x1)*y/2/np.pi)
    return res

th = get_theta(z, z0=z0, z1=z1)
inds = (50,50)
fig, axs = plt.subplots(1,2, figsize=(10,5))
axs[0].imshow((th - th[inds] + np.pi) % (2*np.pi) - np.pi)
N = 2
psi0 = get_psi(z, Nx=N, Ny=N)
axs[1].imshow(np.angle(psi0/psi0[inds]))
```

```{code-cell}
z[inds]
```

```{code-cell}
N = 20
inds = (50,50)
psi0 = get_psi(z, Nx=N, Ny=N)
plt.plot((f[:, 50]/f[inds]).real)
plt.plot((psi0[:, 50]/psi0[inds]).real)
```

```{code-cell}

```

```{code-cell}
inds = (25,25)
th = get_theta(z, z0=z0, z1=z1)
f = np.exp(1j*th)
plt.imshow(abs(f/f[inds] -  psi0/psi0[inds]))
plt.colorbar()
```

```{code-cell}
1.*(abs(z) < 1) - (abs(z) < 2)
```

```{code-cell}

```