{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  $N$ interacting masses in $\\cos$ potential\n",
    "Additional requirements for this example: `scipy`\n",
    "\n",
    "### Sine Gordon equation\n",
    "\n",
    "Let $x$ be an array with position of $N$ masses in  $\\cos$ potential  linearly coupled to nearest neighbours.  Then the Newton's equations read:\n",
    "\n",
    "$$ \\begin{cases}\n",
    "\\displaystyle\\frac{dx}{dt} = v \\\\\n",
    "\\displaystyle\\frac{dv}{dt} = - M x + \\sin(x),\\;\\;\\;(1)\n",
    "\\end{cases} $$\n",
    "\n",
    "where M is a discrete Laplacian. \n",
    "\n",
    "On the other hand this is an aproximation to a sine-Gordon equation: \n",
    "\n",
    "$$ x_{tt}- x_{\\xi\\xi} + \\sin x= 0 $$\n",
    "\n",
    "which has the following soliton solutions:\n",
    "\n",
    "$$x_\\text{soliton}(\\xi, t) := 4 \\arctan \\exp\\left\\{ \\frac{ \\xi - v t }{\\sqrt{1 - v^2}}\\right\\}$$\n",
    "\n",
    "Below we solve numerically the system (1) and animate results in 3D.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 153\n",
    "M = np.diag((N-1)*[ 1.0],-1)+np.diag((N-1)*[ 1.0],1)+np.diag(N*[-2.0],0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.zeros(2*N)\n",
    "x = y0[:N] # Fist N variables are positions\n",
    "v = y0[N:]\n",
    "ksi = np.linspace(-5,25,N, dtype=np.float32)\n",
    "h = np.diff(ksi)[0]\n",
    "\n",
    "instanton = lambda x,v,t: 4 * np.arctan(np.exp( (x-v*t)/np.sqrt(1-v**2) ))\n",
    "v1,v2 = 0.4,0.05\n",
    "x[:] = instanton(ksi,v1,0)\n",
    "v[:] = -v1/h*np.gradient( instanton(ksi,v1,0) )\n",
    "\n",
    "x[:] += instanton(ksi,v2,50)\n",
    "v[:] += -v2/h*np.gradient( instanton(ksi,v2,50) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 982 ms, sys: 197 ms, total: 1.18 s\n",
      "Wall time: 315 ms\n"
     ]
    }
   ],
   "source": [
    "def lhs(y_,t):\n",
    "    y = y_.copy()\n",
    "    x = y_[:N]\n",
    "    v = y_[N:]\n",
    "    y[:N] = v \n",
    "    y[N:] = -np.sin(x) + 1/h**2* np.dot(M,x)  \n",
    "    y[0] = 0\n",
    "    y[N-1] = 0 \n",
    "    return y\n",
    "ts = np.linspace(0,140,50)  \n",
    "%time xt = odeint(lhs,y0, ts).astype(np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ab5b9e228e0e42a38ef522ad8b5bd7fc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import k3d\n",
    "import numpy as np\n",
    "import time \n",
    "plot = k3d.plot()\n",
    "r = .6\n",
    "origins = np.vstack([ksi,np.zeros(N),np.zeros(N)]).T.copy().astype(np.float32)\n",
    "vectors = np.vstack([np.zeros(N),r*np.sin(xt[0,:N]),r*np.cos(xt[0,:N])]  ).T.astype(np.float32)\n",
    "\n",
    "vector_plot = k3d.vectors(origins, vectors)\n",
    "line_plot = k3d.line(vectors + origins,color=0xff0000)\n",
    "\n",
    "plot += vector_plot\n",
    "plot += line_plot\n",
    "plot.display()\n",
    "\n",
    "def update_plot(xx):\n",
    "    vectors = np.vstack([np.zeros_like(xx),r*np.sin(xx),r*np.cos(xx)]  ).T\n",
    "    vector_plot.vectors = vectors.copy()\n",
    "    line_plot.vertices = vectors+origins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot.camera_auto_fit = False\n",
    "plot.grid_auto_fit = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 137 ms, sys: 15.1 ms, total: 152 ms\n",
      "Wall time: 2.73 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "import time\n",
    "plot.auto_rendering = False\n",
    "for xx in xt[:,:N]:\n",
    "    update_plot(xx)\n",
    "    plot.render()\n",
    "    time.sleep(0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5104a59924924abfb7903b69750c08d1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='ith', max=49), Output()), _dom_classes=('widget-interact…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import widgets,interact\n",
    "\n",
    "@interact(ith = widgets.IntSlider(min=0,max=(ts.size-1)))\n",
    "def draw(ith):\n",
    "    update_plot(xt[ith,:N])\n",
    "    plot.render()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  },
  "nbTranslate": {
   "displayLangs": [
    "en",
    "pl"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "pl",
   "targetLang": "en",
   "useGoogleTranslate": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}