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

```{code-cell}
%pylab inline --no-import-all
from IPython.display import display, clear_output
from mmfutils.contexts import NoInterrupt
from pytimeode.evolvers import EvolverABM, EvolverSplit
from gpe.minimize import MinimizeState

import gpe.bec

reload(gpe.bec)
from gpe import utils

from gpe.bec import u


class State(gpe.bec.State):
    def __init__(self, Nx=4800 / 4, Lx=320 * u.micron / 1.5):
        self.U0 = 39 * u.nK
        self.Us = -6 * u.nK
        self.mu = 8.0 * u.nK
        self.v_s = 0.21 * u.mm / u.s
        self.x0 = -80 * u.micron
        self.w0 = 5 * u.micron
        self.lam = 812 * u.nm
        gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(Lx,))

    @property
    def x_s(self):
        return self.x0 + self.t * self.v_s

    def get_Vext(self):
        x = self.xyz[0]
        x0 = np.pi * self.w0**2 / self.lam
        w = self.w0 * np.sqrt(1 + (x / x0) ** 2)
        d_x = 1 * u.micron
        V_s = self.Us * np.vectorize(utils.step)(self.x_s - x - d_x / 2, d_x)
        return self.U0 * (1 - (self.w0 / w) ** 2) + V_s - self.mu


s0 = State()
m = MinimizeState(s0, fix_N=False)
assert m.check()
s = m.minimize()
s.plot()
# plt.plot(s.xyz[0], s.get_Vext())
```

```{code-cell}
e = EvolverABM(s, dt=0.3 * s.t_scale)
# e = EvolverSplit(s, dt=5*s.t_scale)
fig = None
with NoInterrupt() as interupted:
    while not interupted:
        e.evolve(500)
        plt.clf()
        e.y.plot()
        plt.axvline(e.y.x_s)
        display(plt.gcf())
        clear_output(wait=True)
```

# Dumb Holes

```{code-cell}
%pylab inline --no-import-all
from IPython.display import display, clear_output
from mmfutils.contexts import NoInterrupt
from pytimeode.evolvers import EvolverABM, EvolverSplit
from gpe.minimize import MinimizeState

import gpe.soc_soliton

reload(gpe.soc_soliton)
from gpe import utils

from gpe.soc_soliton import u


experiment = gpe.soc_soliton.Experiment()


class State(gpe.soc_soliton.State2tube):
    def __init__(self, experiment=experiment, Nx=4800 / 2, Lx=320 * u.micron):
        self.U0 = 39 * u.nK
        self.Us = -6 * u.nK
        self.mu = 8.0 * u.nK
        self.v_s = 0.21 * u.mm / u.s
        self.x0 = -100 * u.micron
        self.w0 = 5 * u.micron
        self.lam = 812 * u.nm
        gpe.soc_soliton.State2tube.__init__(
            self, experiment=experiment, Nxyz=(Nx,), Lxyz=(Lx,)
        )

    @property
    def x_s(self):
        return self.x0 + self.t * self.v_s

    def get_Vtrap(self):
        x = self.xyz[0]
        x0 = np.pi * self.w0**2 / self.lam
        w = self.w0 * np.sqrt(1 + (x / x0) ** 2)
        d_x = 1 * u.micron
        V_s = self.Us * np.vectorize(utils.step)(self.x_s - x - d_x / 2, d_x)
        V_a = V_b = self.U0 * (1 - (self.w0 / w) ** 2) + V_s - self.mu
        return V_a, V_b

    def plot(self, fig=None):  # pragma: nocover
        from matplotlib import pyplot as plt
        from matplotlib.gridspec import GridSpec

        # from mmfutils.plot import imcontourf

        if fig is None:
            fig = plt.figure(figsize=(15, 3))

        if self.dim == 1:
            gs = GridSpec(3, 1, width_ratios=[5], hspace=0.0)
            x = self.xyz[0] / u.micron
            na, nb = self.get_density()
            t_ = self.t / self.t_unit
            E_R = self.experiment.E_R

            plt.subplot(gs[:, 0])
            plt.plot(x, na, "b:")
            plt.plot(x, nb, "g:")
            plt.plot(x, na + nb, "k-")
            plt.xlabel("x [micron]")
            plt.tight_layout()
        else:
            raise NotImplementedError("Only 1D plotting supported")

        E = self.get_energy()
        Na, Nb = self.get_Ns()
        N = Na + Nb
        plt.suptitle(
            "t={}t0, Ns={:.4f}+{:.4f}={:.4f}, E={:.4f}".format(
                self.t / self.t0, Na, Nb, N, E
            )
        )
        return fig


s0 = State()
s0.plot()
m = MinimizeState(s0, fix_N=False)
assert m.check()
s = m.minimize()
s.plot()
# plt.plot(s.xyz[0], s.get_Vext())
```