---
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}
:init_cell: true

import mmf_setup;mmf_setup.nbinit()
%pylab inline --no-import-all
```

# Utils

+++

This notebook describes some of the utilities provided here.

+++

## Periodic Potential Extension

+++

A common case is that of a large trap along the $x$ direction where we are interested in the behavior near the center.  To enable this case, we provide a flag `cells_x` which allows the user to specify a small number of $k_r$ cells near the center of the trap to simulate.  To avoid numerical issues at the boundaries, we need a transformation so that the potential can be applied in such a way as to be $C_{\infty}$ periodic over this new interval.  There are several options:

+++

### Harmonic Potentials

If the potential is Harmonic $V(x) \propto x^2$, then we can replace it with a periodic function that matches near the center.  This is done by the function [`gpe.utils.x2_2()`](../gpe/utils.py:55) which is a periodic function $f(x)$ with period $2\pi$ approximating $x^2/2$ at the center:

$$
  V(x) = \frac{m\omega^2}{2}\frac{x^2}{2}, \qquad
  V_{\mathrm{periodic}}(x) = \frac{m\omega^2}{2 k^2}f_{\mathrm{order}}(kx), \qquad
  k = \frac{2\pi}{L}
$$

\begin{align}
    f_0(x) &= 1 - \cos x = \frac{x^2}{2} + \order(x^4),\\
    f_2(x) &= \frac{1}{6}\cos^2 x - \frac{4}{3}\cos x + \frac{7}{6} = \frac{x^2}{2} + \order(x^6),\\
    f_4(x) &= \frac{-2}{45}\cos^3 x + \frac{3}{10}\cos^2 x - \frac{22}{15}\cos x + \frac{109}{90} = \frac{x^2}{2} + \order(x^8),\\
    f_6(x) &= \frac{1}{70}\cos^4 x - \frac{32}{315}\cos^3 x + \frac{27}{70}\cos^2 x - \frac{32}{21}\cos x + \frac{386}{315} = \frac{x^2}{2} + \order(x^{10}).
\end{align}

```{code-cell}
from gpe.utils import x2_2

x = np.linspace(-np.pi, np.pi, 100)
for order in [0,2,4,6]:
    plt.plot(x, x2_2(x, order=order), label="order={}".format(order))
plt.plot(x, x**2/2, ':', label="$x^2/2$", scaley=False)
plt.legend(); plt.xlabel('$x$'); plt.ylabel('x2_x(x, order)')
    
```

### Abscissa Transform

+++

Another option is to transform the abscissa so as to render the potential periodic.  If $V(x) = V(-x)$, then the following transform works quite nicely, rendering the potential periodic over the interval $x\in [-1,1]$:

$$
  V_{\mathrm{periodic}}(x) = V(\tilde{x}), \qquad
  \tilde{x} = \left[\frac{1}{a}\tanh\left(\frac{2 a}{\pi}\tan\frac{\pi x^p}{2}\right)\right]^{1/p}
$$

```{code-cell}
import gpe.utils;reload(gpe.utils)
from gpe.utils import x_periodic

def V(x):
    return x**2/2

def Vp(x, x0, p=3):
    xt = x_periodic(x, x0=x0, p=p)
    return V(xt)

x = np.linspace(-1, 1, 100)
for x0 in [0.8, 0.9, 1.0]:
    plt.plot(x, Vp(x, x0=x0, p=3), label="x0={}".format(x0))
plt.plot(x, V(x), ':', label="$V(x)$", scaley=False)
plt.legend(); plt.xlabel('$x$'); plt.ylabel('V(x, p)')
```

Here is a comparison of the two methods with the default values (chosen to match roughly):

```{code-cell}
from gpe.utils import x2_2, x_periodic


w = 1.5
m = 2.3
L = 3.0

def V(x):
    return m*(w*x)**2/2

def Vp1(x):
    k = 2*np.pi/L
    return m*w**2/k**2 * x2_2(k*x)

def Vp2(x):
    xt = L/2*x_periodic(2*x/L)
    return V(xt)

x = np.linspace(-L/2, L/2, 100)

# This is the transform which shifts the potential to
# be centered at x0
x0 = 1.0
x_ = (x - x0 - x.min()) % L + x.min()
plt.plot(x, V(x_))
plt.plot(x, Vp1(x_), label="Cosine Replacement")
plt.plot(x, Vp2(x_), label="Periodic Extension")
plt.legend()
```

```{code-cell}
        Lx = self._get_Lx(state)
        xyz[0] = (x - self.get('x0', t_=t_) - x.min()) % Lx + x.min()
        return xyz
```