{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Paraboloidal pendulum\n",
    "\n",
    "## System definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t')\n",
    "var('l g')\n",
    "xy_names = [('x','x'),('y','y'),('z','z')]\n",
    "load('cas_utils.sage')\n",
    "to_fun, to_var = make_symbols(xy_names)\n",
    "xy = [vars()[v] for v,lv in xy_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dAlemb = (x.subs(to_fun).diff(t,2))*dx + \\\n",
    "         (y.subs(to_fun).diff(t,2))*dy + \\\n",
    "         (z.subs(to_fun).diff(t,2)+g)*dz  \n",
    "dAlemb = dAlemb.subs(to_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(dAlemb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 1/2*(x^2+y^2)-z\n",
    "dxy = [vars()['d'+repr(zm)] for zm in xy]\n",
    "\n",
    "constr =sum([dzm*f.diff(zm) for zm,dzm in zip(xy,dxy)])\n",
    "showmath(constr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq1=(dAlemb.subs(constr.solve(dz)[0])*x).expand().coefficient(dx).subs(to_var)\n",
    "eq2=(dAlemb.subs(constr.solve(dz)[0])*x).expand().coefficient(dy).subs(to_var)\n",
    "\n",
    "showmath([eq1,eq2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve([f.subs(to_fun).diff(2).subs(to_var),eq1,eq2],[xdd,ydd,zdd])[0][:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode = [xd,yd] + [s_.rhs() for s_ in sol]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode = map(lambda x:x.subs({l:1,g:1}),ode)\n",
    "showmath(ode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = srange(0,237,.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_prime = (1/2*(x^2+y^2)).subs(to_fun).diff(t).subs(to_var)\n",
    "z_prime.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin = (1/2*(xd^2+yd^2+z_prime^2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numsol = desolve_odeint(ode,[1.0,0,0.,.5],times,[x,y,xd,yd])\n",
    "p = line ( zip(numsol[:,0],numsol[:,1]) )\n",
    "Epot = 1/2*(numsol[:,0]^2+numsol[:,1]^2)\n",
    "p += circle( (0,0), sqrt(np.max(2*Epot)),color='red' )\n",
    "p += circle( (0,0), sqrt(np.min(2*Epot)),color='green' )\n",
    "p.show(aspect_ratio=1,axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin_num = [Ekin.subs({x:d[0],y:d[1],xd:d[2],yd:d[3]}) for d in numsol]\n",
    "Epot_num  = 1/2*(numsol[:,0]^2+numsol[:,1]^2) # g=1!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin_num + Epot_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Angular momentum\n",
    "\n",
    "We can calculate symbolically angular momentum and chech if its $z$-component is conserved in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x y', domain='real')\n",
    "\n",
    "#e_r = vector([x,y,0])\n",
    "#v = vector([xd,yd,0])\n",
    "\n",
    "e_r = vector([x,y,1/2*(x^2+y^2)])\n",
    "v = vector([xd,yd,z_prime])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = e_r.cross_product(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "showmath(p[0].full_simplify())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin_num = [Ekin.subs({x:d[0],y:d[1],xd:d[2],yd:d[3]}) for d in numsol]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_num = [(p[2]).subs({x:d[0],y:d[1],xd:d[2],yd:d[3]}) for d in numsol]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_num[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin.subs( (p[2]*x*y).expand().solve(xd)[0]*x).expand().show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ekin.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_num[:4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact \n",
    "def plot2d_traj(v0=slider(0,3,0.001)):\n",
    "    numsol = desolve_odeint(ode,[1,0,0,v0],times,[x,y,xd,yd])\n",
    "    p = line( zip(numsol[:,0],numsol[:,1]), aspect_ratio=1)\n",
    "    p += circle( (0,0),1,color='red')\n",
    "    p.show(axes=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot2d_traj(v0=0.23),plot2d_traj(v0=0.63),plot2d_traj(v0=1.63),"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3d = line3d( zip(numsol[:,0],numsol[:,1],(1/2*(numsol[:,0]**2+numsol[:,1]**2))),thickness=2,color='red')\n",
    "p3d += plot3d(1/2*(x^2+y^2),(x,-1.2,1.2),(y,-1.2,1.2),opacity=0.8) \n",
    "p3d.show(viewer='tachyon')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\newpage"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.7",
   "language": "",
   "name": "sagemath"
  },
  "language": "python",
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "nbTranslate": {
   "displayLangs": [
    "en",
    "pl"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "pl",
   "targetLang": "en",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}