---
execution:
  timeout: 60
jupytext:
  cell_metadata_filter: -all
  formats: md:myst,ipynb
  main_language: python
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.7
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:tags: [hide-cell]

import mmf_setup; mmf_setup.nbinit()
import os, sys
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None

from matplotlib.animation import FuncAnimation
from gpe.contexts import FPS
```

(sec:Viscosity1)=
# Viscosity Notes 1

These are some supporting notes for the document {ref}`sec:EffectiveViscosity`.

## Numerical Potentials

In {ref}`sec:BarrierFlow`, we want some potentials that have a barrier height $V_0$, a
width $\sigma$, and a drop $\delta$.  This is easy to do with a splines:

```{code-cell}
from scipy.interpolate import InterpolatedUnivariateSpline

def V_linear(x, V0=0.5, sigma=2.0, dV=0.2, d=0, k=1):
    """Return the dth derivative of a linear potential barrier.
    
    Arguments
    ---------
    x : array-like
        Abscissa on which to compute the potential.
    V0 : float
        Height of barrier.
    sigma : float
        Width of barrier.
    dV : float
        Potential difference from right to left of the barrier.
    d : [0, 1]
        Derivative order.
    k : int
        Order of spline.
    """
    xp, fp = [0, sigma / 2, sigma], [0, V0, dV]
    if k > 1:
        xp = sigma * np.linspace(0, 1, 7)
        fp = [0] * 3 + [V0] + [dV] * 3
        xp = sigma * np.array([0, 0.001, 0.5, 0.999, 1])
        fp = [0] * 2 + [V0] + [dV] * 2
    V = InterpolatedUnivariateSpline(xp, fp, k=k, ext="const")
    while d > 0:
        V = V.derivative()
        d -= 1
    return V(x)
    

V0 = 1.0
sigma = 1.0
dV = 0.5

x = sigma * np.linspace(-0.5, 1.5, 101)

fig, ax = plt.subplots(figsize=(4, 3))
ax1 = ax.twinx()

for k in [1, 2, 3]:
    V = V_linear(x, sigma=sigma, V0=V0, dV=dV, k=k)
    dV_dx = V_linear(x, sigma=sigma, V0=V0, dV=dV, k=k, d=1)
    ax.plot(x / sigma,  V / V0, label=f"{k=}")
    ax1.plot(x / sigma,  dV_dx / (V0/sigma), ":")

ax.legend()
ax.set(xlabel=r"$x/\sigma$", ylabel="$V/V_0$",
       xticks=[-0.5, 0, 0.5, 1.0, 1.5])
ax1.set(ylabel=r"$V'(x)$ [$V_0/\sigma$]");
```

We can also do this with a gaussian and the error function if we want something
analytic.

```{code-cell}
from scipy.special import erf

def V_gaussian(x, V0=0.5, sigma=2.0, dV=0.2, d=0, V0_supp=10.0):
    """Return the dth derivative of a linear potential barrier.
    
    Arguments
    ---------
    x : array-like
        Abscissa on which to compute the potential.
    V0 : float
        Height of barrier.
    sigma : float
        Width of barrier.
    dV : float
        Potential difference from right to left of the barrier.
    d : [0, 1]
        Derivative order.
    V0_supp : float
        The potential is shifted so that `V(0) = exp(-V_0supp)`.
    """
    s = sigma / 2
    f2 = V0_supp
    f = np.sqrt(f2) 
    V0_ = V0 * np.exp(- f2 * (x / s - 1) ** 2)
    V1_ = dV * (1 + erf( f * (x / s - 1)))/2
    if d == 0:
        return V0_ + V1_
    elif d == 1:
        dV0_ = -V0 * 2 * f2 / s
        dV1_ = dV / np.sqrt(np.pi) * f / s
        return  (dV0_ * (x / s - 1) + dV1_) * np.exp(- f2 * (x / s - 1) ** 2)
    
V0 = 1.0
sigma = 1.0
dV = 0.5

x = sigma * np.linspace(-0.5, 1.5, 301)

V = V_gaussian(x, sigma=sigma, V0=V0, dV=dV)
dV_dx = V_gaussian(x, sigma=sigma, V0=V0, dV=dV, d=1)
assert np.allclose(np.gradient(V, x), dV_dx, rtol=0.02)

fig, ax = plt.subplots(figsize=(4, 3))
ax1 = ax.twinx()
ax.plot(x / sigma,  V / V0)
ax1.plot(x / sigma,  dV_dx / (V0/sigma), ":")

ax.set(xlabel=r"$x/\sigma$", ylabel="$V/V_0$",
       xticks=[-0.5, 0, 0.5, 1.0, 1.5])
ax1.set(ylabel=r"$V'(x)$ [$V_0/\sigma$]");
```