ball-bouncing-uneven-street-HC-method.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sm\n",
    "import sympy.physics.mechanics as me\n",
    "import time\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "from scipy.optimize import fsolve, minimize\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mp\n",
    "from matplotlib import patches\n",
    "%matplotlib inline\n",
    "from IPython.display import HTML\n",
    "mp.rcParams['animation.embed_limit'] = 2**128\n",
    "from matplotlib import animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Needed to only plot a few *hysteresis curves*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class Rausspringen(Exception):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An elastic homogenious ball of radius *r* and mass *m* is thrown (or dropped) on an uneven street.\\\n",
    "An observer (a particle) of mass $m_o$ may be put anywhere inside the ball.\n",
    "\n",
    "**Variables**\n",
    "\n",
    "- $N$: inertial frame\n",
    "- $A$: frame fixed to the ball\n",
    "\n",
    "- $Dmc$: center of the ball, a *RigidBody*\n",
    "- $P_o$: the observer, a *Particle*\n",
    "- $P_{st}$: point where the ball hits the street.\n",
    "\n",
    "- $q_1, q_2, q_3$: angle of the ball w.r.t. N\n",
    "- $u_1, u_2, u_3$: angular velocities of the ball\n",
    "- $x_b, y_b, z_b$: location of the center of the ball, relative to the inertial frame N\n",
    "- $ux_b, uy_b, uz_b$: the speeds\n",
    "\n",
    "\n",
    "- $c_{tau}$: the experimental constant needed for Hunt-Crossley\n",
    "- $m_u$: the friction constant\n",
    "- $ny_b, ny_s$: Poisson's ratio of ball and street respectively.\n",
    "- $E_b, E_s$: Young's modulus of ball and street respectively\n",
    "- $rhodt_{street}$: speed right before the impact.\n",
    "\n",
    "\n",
    "- $x_{st}, z_{st}$: impact point, where the ball hits the street.\n",
    "\n",
    "- $m, m_o, i_{XX}, i_{YY}, i_{ZZ}$: mass of ball, mass of observer, moments of inertia of the ball\n",
    "- $amplitude, frequenz$: parameters to describe the bumpy road\n",
    "- $\\alpha, \\beta, \\gamma$: describe the location of the observer relative to the center of the ball \n",
    "\n",
    "NOTE: This seems very 'tricky' as far as the numerical integration is concerned. I assume, because the systen is extremely stiff during collisions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time.time()\n",
    "\n",
    "N, A = sm.symbols('N, A', cls = me.ReferenceFrame)\n",
    "P0, Dmc, Po, Pst = sm.symbols('P0, Dmc, Po, Pst', cls = me.Point)\n",
    "\n",
    "q1, q2, q3, u1, u2, u3 = me.dynamicsymbols('q1, q2, q3, u1, u2, u3')\n",
    "xst, zst = sm.symbols('xst, zst')\n",
    "xb, yb, zb, uxb, uyb, uzb = me.dynamicsymbols('xb, yb, zb, uxb, uyb, uzb')\n",
    "\n",
    "m, mo, g, r, k, amplitude, frequenz, alpha, beta, gamma, iXX, iYY, iZZ, t = sm.symbols('m, mo, g, r,' \n",
    "        'k, amplitude, frequenz, alpha, beta, gamma, iXX, iYY, iZZ, t')\n",
    "\n",
    "nyb, nys, Eb, Es, rhodtstreet, ctau, mu = sm.symbols('nyb, nys, Eb, Es, rhodtstreet, ctau, mu')\n",
    "\n",
    "A.orient_body_fixed(N, (q1, q2, q3), '123')\n",
    "rot = A.ang_vel_in(N)\n",
    "A.set_ang_vel(N, u1*N.x + u2*N.y + u3*N.z)\n",
    "rot1 = A.ang_vel_in(N)\n",
    "\n",
    "P0.set_vel(N, 0.)\n",
    "\n",
    "Dmc.set_pos(P0, xb*N.x + yb*N.y + zb*N.z)\n",
    "Dmc.set_vel(N, uxb*N.x + uyb*N.y + uzb*N.z)\n",
    "\n",
    "Po.set_pos(Dmc, r * (alpha*A.x + beta*A.y + gamma*A.z))\n",
    "Po.v2pt_theory(Dmc, N, A);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Here the street is modeled*.\\\n",
    "Basically it is a parabola open to the top (that is the positive Y direction), with superimposed sine waves, of ever larger frequency, but smaller amplitude.\\\n",
    "The integer *rumpel*, larger than zero, determines how many sine waves will be used.\\\n",
    "At the bottom, I define $P_{st}$, the point of impact of the ball with the street. How I find it is described further down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#============================================\n",
    "rumpel = 3                     \n",
    "#============================================\n",
    "if isinstance(rumpel, int) == False or rumpel < 1:\n",
    "    raise Exception('rumpel must be an integer larger than zero')\n",
    "    \n",
    "def gesamt(x, z, amplitude, frequenz):\n",
    "    strasse = sum([amplitude/j * (sm.sin(j*frequenz * x) + 2.5 * sm.sin(j*frequenz * z)) for j in range(1, rumpel)])\n",
    "    strassen_form = (frequenz/2. * x)**2  + (frequenz/2. * z)**2\n",
    "    return strassen_form + strasse \n",
    "\n",
    "# needed for the lambdification, as a function cannot be lambdified - at least I do not know how.\n",
    "gesamt1 = gesamt(xb, zb, amplitude, frequenz) \n",
    "\n",
    "Pst.set_pos(P0, xst*N.x + gesamt(xst, zst, amplitude, frequenz)*N.y + zst*N.z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the minimal osculating circle (Schmiegekreis) of the street. I found the formula for the one dimensional case in the internet.\\\n",
    "I simply take the smaller ones of the X - direction, Z - direction, respectively.\\\n",
    "There may the a formula for this 3D case, I did not bother to find it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_max_xst = (sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(xst))**2\n",
    "           )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(xst, 2)\n",
    "r_max_zst = (sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(zst))**2 \n",
    "           )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(zst, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Point of impact**\\\n",
    "Here I get the $x_{st}, z_{st}$ coordinates of the point where the balls hits the street.\\\n",
    "The idea is this:\\\n",
    "$vector_1$ is normal to the street at the point $P_{st}$ = ($x_{st}$, gesamt($x_{st}, z_{st}$, *amplitude, frequenz*), $z_{st}$). It points 'inwards'. I got this formula from the internet.\\\n",
    "$vector_2$ is the vector pointing from $P_{st}$ to the center of the ball, *Dmc*\\\n",
    "These two vectors must be parallel at all times, that is $\\hat{vector_1}$ = $\\hat{vector_2}$, with $\\hat a$ meaning $\\dfrac{a}{|a|}$. This means, that $P_{st}$ would be the contact point of the ball with the street, if its radius *r* was equal to |$vector_2|$.\\\n",
    "This gives two equations for $x_{st}, z_{st}$ to be solved numerically later, to get the correct starting values of $x_{st}, z_{st}$, and then again in every step of the numerical integration.\\\n",
    "This idea works at least as long as there is only one solution to the equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vektor1 = -(gesamt(xst, zst, amplitude, frequenz).diff(xst)*N.x - N.y + \n",
    "           gesamt(xst, zst, amplitude, frequenz).diff(zst)*N.z).normalize()\n",
    "\n",
    "vektor2 = Dmc.pos_from(Pst).normalize()\n",
    "\n",
    "auftreff = vektor1 - vektor2\n",
    "\n",
    "auftreffX = me.dot(auftreff, N.x)\n",
    "auftreffZ = me.dot(auftreff, N.z)\n",
    "\n",
    "loesung = sm.Matrix([auftreffX, auftreffZ])\n",
    "print('loesung DS', me.find_dynamicsymbols(loesung))\n",
    "print('loesung free symbols', loesung.free_symbols)\n",
    "print('loesung has {} operations'.format(sum([loesung[i].count_ops(visual=False) for i in range(2)])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Function which calculates the forces of an impact between the ball and the street**\\\n",
    "1.\\\n",
    "The impact force is on the line normal to the street, going through the center of the ball, $P_1$. I use Hunt Crossley's method to calculate it\\\n",
    "It's direction points from the street to the ball.\\\n",
    "2.\\\n",
    "The friction force acting on the contact point $CP_0$ is proportional to the component of the speed of $CP_0$ in the plane tanget to the street, the impact force and a friction factor $m_{uW}$.\\\n",
    "It is equivalent to a force through $P_1$ and to a torque acting on $A_1$, the ball fixed frame.\n",
    "\n",
    "\n",
    "**Note about the force during the collisions**\n",
    "\n",
    "**Hunt Crossley's method**\n",
    " \n",
    "My reference is this article, given to me by JM\\\n",
    "https://www.sciencedirect.com/science/article/pii/S0094114X23000782 \\\n",
    "\n",
    " \n",
    "This is with dissipation during the collision, the general force is given in (63) as\\\n",
    "$f_n = k_0 \\cdot \\rho + \\chi \\cdot \\dot \\rho$, with $k_0$ as above, $\\rho$ the penetration, and $\\dot\\rho$ the speed of the penetration.\\\n",
    "In the article it is stated, that $n = \\frac{3}{2}$ is a good choice, it is derived in Hertz' approach. Of course, $\\rho, \\dot\\rho$ must be the signed magnitudes of the respective vectors.\n",
    "\n",
    "A more realistic force is given in (64) as:\\\n",
    "$f_n = k_0 \\cdot \\rho^n + \\chi \\cdot \\rho^n\\cdot \\dot \\rho$, as this avoids discontinuity at the moment of impact.\n",
    "\n",
    "**Hunt and Crossley** give this value for $\\chi$, see table 1:\n",
    "\n",
    "$\\chi = \\dfrac{3}{2} \\cdot(1 - c_\\tau) \\cdot \\dfrac{k_0}{\\dot \\rho^{(-)}}$, \n",
    "where $c_\\tau = \\dfrac{v_1^{(+)} - v_2^{(+)}}{v_1^{(-)} - v_2^{(-)}}$, where $v_i^{(-)}, v_i^{(+)}$ are the speeds of $body_i$, before and after the collosion, see (45), $\\dot\\rho^{(-)}$ is the speed right at the time the impact starts. $c_\\tau$ is an experimental factor, apparently around 0.8 for steel.\n",
    "\n",
    "Using (64), this results in their expression for the force:\n",
    "\n",
    "$f_n = k_0 \\cdot \\rho^n \\left[1 + \\dfrac{3}{2} \\cdot(1 - c_\\tau) \\cdot \\dfrac{\\dot\\rho}{\\dot\\rho^{(-)}}\\right]$\n",
    "\n",
    "with $k_0 = \\frac{4}{3\\cdot(\\sigma_1 + \\sigma_2)} \\cdot \\sqrt{\\frac{R_1 \\cdot R_2}{R_1 + R_2}}$, where $\\sigma_i = \\frac{1 - \\nu_i^2}{E_i}$, with $\\nu_i$ = Poisson's ratio, $E_i$ = Young\"s modulus, $R_1, R_2$ the radii of the colliding bodies, $\\rho$ the penetration depth. All is near equations (54) and (61) of this article.\\\n",
    "From the ball's perspective, the street is concave. In the formula I simply take $R_2 < 0$. As the street has different osculating radii, I take the average of the radii in the X / Z directions. As |osculating radius| > |radius of the ball|, the term of the square root should always be positive.\\\n",
    "I am not sure, whether this applies the formula correctly.\n",
    "\n",
    "As per the article, $n = \\frac{3}{2}$ is always to be used.\n",
    "\n",
    "*spring energy* =   $ k_0 \\cdot \\int_{0}^{\\rho} k^{3/2}\\,dk$ = $k_0 \\cdot\\frac{2}{5} \\cdot \\rho^{5/2}$\\\n",
    "I assume, the dissipated energy cannot be given in closed form, at least the article does not give one.\n",
    "\n",
    "*Note*  \n",
    "$c_\\tau = 1.$ gives **Hertz's** solution to the impact problem, also described in the article.\n",
    "\n",
    "\n",
    "**Friction when the ball hits the street**\n",
    "\n",
    "This website:\n",
    "https://math.stackexchange.com/questions/2195047/solve-the-vector-cross-product-equation\n",
    "\n",
    "gives: $b = a \\times x \\rightarrow x = \\dfrac{b \\times a}{|a|^2}$\n",
    "\n",
    "This way, I can easily get the force acting on CP0 ( = contact point of ball with wall), without any further geometric considerations. Of course the friction force has opposite to the speed of CP0\\\n",
    "The friction force on CP0 is equivalent to a force on P1 (center of the ball) and a torque on A1 (ball fixed frame)\n",
    "\n",
    "In the formula for the impact force, *forcec* I use |$rhodtstreet$| and  **-**$rhodt$. The reason is this:\n",
    "- As the diagram below shows, the speed jumps dramatically during impact. Hence when I try to get $\\dot \\rho^{(-)}$ it may be positive or negative.\n",
    "- As per my definition the speed component in the first phase of the impact is negative. Hence during the first phase of the impact, $\\dfrac{\\dot \\rho}{\\dot \\rho^{(-)}} > 0.$, as it should be.\n",
    "\n",
    "*vorzeichen* is introduced to ensure, that the collision force continues to act if the ball falls through the street."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HC_street(N, A1, P1, r, ctau, rhodtstreet):\n",
    "        \n",
    "    vektor = P1.pos_from(Pst)\n",
    "    vorzeichen = sm.sign(me.dot(vektor, N.y)) # positive if the ball is above the street, lse negative.\n",
    "    abstand = vektor.magnitude()\n",
    "    richtung = vektor.normalize()\n",
    "    rho = r - abstand       # positive in the ball has penetrated the wall\n",
    "    CP0 = me.Point('CP0')\n",
    "    CP0.set_pos(P1, -vektor)\n",
    "    vCP0 = CP0.v2pt_theory(Dmc, N, A1)\n",
    "    rhodt = me.dot(vCP0, richtung)             #me.dot(P1.pos_from(P0).diff(t, N), richtung)\n",
    "    rho = sm.Abs(rho) # sm.Max(rho, sm.S(0))\n",
    "\n",
    "#determine R2 of the a.m. formula    \n",
    "    r_xst   = -(sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(xst))**2\n",
    "               )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(xst, 2)\n",
    "    r_zst   = -(sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(zst))**2 \n",
    "               )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(zst, 2)\n",
    "    R2      = 0.5 * (r_xst + r_zst)\n",
    "# determine k0\n",
    "# the sm.Abs(..) should theoretically not be necessary. But seems to avoid issues at times.\n",
    "    wurzel = sm.sqrt(sm.Abs(r * R2 /(r + R2)))\n",
    "    sigma_b = (1. - nyb**2) / Eb\n",
    "    sigma_s = (1. - nys**2) / Es\n",
    "    k0      = 4. / (3.*(sigma_b + sigma_s)) * wurzel\n",
    "\n",
    "# Here I assume, that the penetration is small, otherwise I would not know how to do it.\n",
    "    forcec = k0 * rho**(3/2) * ( 1. + 3./2. * (1 - ctau) * (-rhodt) / sm.Abs(rhodtstreet)) * (richtung *\n",
    "            sm.Heaviside(r - vorzeichen * abstand, 0.))\n",
    "    \n",
    "    vx = vCP0 - (me.dot(vCP0, richtung)) * richtung\n",
    "    friction_force = forcec.magnitude() * mu * (-vx)        \n",
    "    hilfs = CP0.pos_from(P1)\n",
    "    torque = hilfs.cross(friction_force) * sm.Heaviside(r - abstand, 0.)\n",
    "    forcef = 1./me.dot(hilfs, hilfs) * torque.cross(hilfs) * sm.Heaviside(r - abstand, 0.)\n",
    "        \n",
    "    return [forcec, forcef, torque]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Kane's equations**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "start2 = time.time()\n",
    "\n",
    "I = me.inertia(A, iXX, iYY, iZZ)\n",
    "body = me.RigidBody('body', Dmc, A, m, (I, Dmc))\n",
    "Poa  = me.Particle('Poa', Po, mo)\n",
    "BODY = [body, Poa]\n",
    "\n",
    "F1 = [(Dmc, -m*g*N.y), (Po, -mo*g*N.y)]\n",
    "F2 = [(Dmc, HC_street(N, A, Dmc, r, ctau, rhodtstreet)[0] + HC_street(N, A, Dmc, r, ctau, rhodtstreet)[1])]\n",
    "F3 = [(A, HC_street(N, A, Dmc, r, ctau, rhodtstreet)[2])]\n",
    "\n",
    "FL = F1 + F2 + F3\n",
    "\n",
    "q_ind = [xb, yb, zb] + [q1, q2, q3]\n",
    "u_ind = [uxb, uyb, uzb] + [u1, u2, u3]\n",
    "\n",
    "kd = [me.dot(rot - rot1, uv) for uv in A] + [uxb - xb.diff(t), uyb - yb.diff(t), uzb - zb.diff(t)]\n",
    "\n",
    "KM = me.KanesMethod(N, q_ind=q_ind, u_ind=u_ind, kd_eqs=kd)\n",
    "(fr, frstar) = KM.kanes_equations(BODY, FL)\n",
    "\n",
    "MM = KM.mass_matrix_full\n",
    "force = KM.forcing_full\n",
    "\n",
    "print('force DS', me.find_dynamicsymbols(force))\n",
    "print('force free symbols', force.free_symbols)\n",
    "print('force has {} operations'.format(sum([force[i].count_ops(visual=False) \n",
    "                for i in range(len(force))])), '\\n')\n",
    "\n",
    "print('MM DS', me.find_dynamicsymbols(MM))\n",
    "print('MM free symbols', MM.free_symbols)\n",
    "print('MM has {} operations'.format(sum([MM[i, j].count_ops(visual=False) \n",
    "                for i in range(MM.shape[0]) for j in range(MM.shape[1])])), '\\n')\n",
    "\n",
    "print('It took {:.3f} sec to establish Kanes equations, for rumpel = {}'.format(time.time() - start2\n",
    "        , rumpel))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here some *functions*, needed further down are defined:\n",
    " - $abstand_1$: needed to calculate the distance of the center of the ball from the street\n",
    " - $rhodtstreet_1$: needed to calculate the speed of the impact.\n",
    " \n",
    " \n",
    "I like to look at the *energies of the system*. This has often (just about always!) told me, I had made a mistake in the equations of motion.\\\n",
    "So, these functions are defined here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pot_energie = m * g * me.dot(Dmc.pos_from(P0), N.y) + mo * g * me.dot(Po.pos_from(P0), N.y)\n",
    "kin_energie = sum([koerper.kinetic_energy(N) for koerper in BODY])\n",
    "\n",
    "rho  = Dmc.pos_from(Pst).magnitude()\n",
    "rho1 = sm.Max(r - rho, sm.S(0)) # penetration depth, positive if there is penetration.\n",
    "\n",
    "#determine R2 of the a.m. formula    \n",
    "r_xst   = -(sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(xst))**2\n",
    "               )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(xst, 2)\n",
    "r_zst   = -(sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(zst))**2 \n",
    "               )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(zst, 2)\n",
    "R2      = 0.5 * (r_xst + r_zst)\n",
    "\n",
    "# determine k0    \n",
    "wurzel = sm.sqrt(sm.Abs(r * R2 /(r + R2))) # theoretically the sm.Abs(..) should not be needed.\n",
    "sigma_b = (1. - nyb**2) / Eb\n",
    "sigma_s = (1. - nys**2) / Es\n",
    "k0      = 4. / (3.*(sigma_b + sigma_s)) * wurzel\n",
    "spring_energie = 2./5. * k0 * sm.Heaviside(r - rho, 0.) * rho1**(5/2)\n",
    "\n",
    "vektor = Dmc.pos_from(Pst)\n",
    "abstand1 = vektor.magnitude()\n",
    "\n",
    "richtung = vektor.normalize()\n",
    "CP0 = me.Point('CP0')\n",
    "CP0.set_pos(Dmc, -vektor)\n",
    "vCP0 = CP0.v2pt_theory(Dmc, N, A)\n",
    "rhodtstreet1 = me.dot(vCP0, richtung)\n",
    "print('vCP0 DS', me.find_dynamicsymbols(vCP0, reference_frame=N))\n",
    "print('vCP0 FS', vCP0.free_symbols(reference_frame=N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the *sympy functions* to *numpy functions* so numerical evaluations may be done. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "qL   = q_ind + u_ind\n",
    "pL   = [xst, zst] + [m, mo, g, r, amplitude, frequenz, alpha, beta, gamma, iXX, iYY, iZZ\n",
    "      ] + [nyb, nys, Eb, Es, ctau, mu] + [rhodtstreet]\n",
    "pL1  = [xst, zst] + [m, mo, g, r, amplitude, frequenz, alpha, beta, gamma, iXX, iYY, iZZ\n",
    "      ] + [nyb, nys, Eb, Es, ctau, mu]\n",
    "pL11 = [m, mo, g, r, amplitude, frequenz, alpha, beta, gamma, iXX, iYY, iZZ\n",
    "      ] + [nyb, nys, Eb, Es, ctau, mu]\n",
    "    \n",
    "MM_lam = sm.lambdify(qL + pL, MM, cse=True)\n",
    "force_lam = sm.lambdify(qL + pL, force, cse=True)\n",
    "\n",
    "gesamt_lam = sm.lambdify([xb, zb, amplitude, frequenz], gesamt1, cse=True)\n",
    "\n",
    "loesung_lam = sm.lambdify([xst, zst] + [xb, yb, zb] + [amplitude, frequenz], loesung, cse=True)\n",
    "\n",
    "r_max_lam = sm.lambdify([xst, zst] + pL11, [r_max_xst, r_max_zst], cse=True)\n",
    "\n",
    "pot_lam    = sm.lambdify(qL + pL, pot_energie, cse=True)\n",
    "kin_lam    = sm.lambdify(qL + pL, kin_energie, cse=True)\n",
    "spring_lam = sm.lambdify(qL + pL, spring_energie, cse=True)\n",
    "\n",
    "abstand1_lam = sm.lambdify(qL + pL, abstand1, cse=True)\n",
    "\n",
    "rhodtstreet1_lam = sm.lambdify(qL + pL1, rhodtstreet1, cse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *parameters* and the *initial conditions* are set here. For their meaning, see above.\\\n",
    "It makes sense to use names close to the names used to set up Kane's equations, but they should **not** be the same ones. The sm.symbols and the me.dynamic symbols would be overwritten - with unintended consequences.\n",
    "\n",
    "For reasons I do not understand, the *no energy loss* case, that is $c_{tau} \\approx 1.$ and $m_u \\approx 0.$ gives numerical problems. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#==============================================================\n",
    "# Set parameters and initial conditions.\n",
    "xb1, yb1, zb1 = 0., 15., 0.\n",
    "uxb1, uyb1, uzb1 = -1., 1., 1.\n",
    "\n",
    "q11, q21, q31 = 0., 0., 0.\n",
    "u11, u21, u31 = 0., 0., 0.\n",
    "\n",
    "m1 = 1.\n",
    "mo1 = 0.1\n",
    "r1 = 1.\n",
    "\n",
    "nyb1  = 0.5   # Poisson\"s ratio for rubber, from the internet\n",
    "nys1  = 0.15  # Poisson's ration for concrete, from the internet\n",
    "Eb1   = 3.e8  # Young's modulus for rubber is really about 3.e9\n",
    "Es1   = 2.e10 # Young's modulus for steel is really about 2.e11\n",
    "\n",
    "ctau1 = 0.85    # experimental constant\n",
    "mu1   = 0.1    # coefficient of friction\n",
    "\n",
    "alpha1, beta1, gamma1 = 0.7, 0., 0.7\n",
    "amplitude1 = 1.\n",
    "frequenz1 = .15\n",
    "\n",
    "intervall = 7.5\n",
    "\n",
    "#===============================================================\n",
    "if alpha1**2 + beta1**2 + gamma1**2 >= 1.:\n",
    "    raise Exception('Observer is outside of the ball')\n",
    "    \n",
    "schritte     = int(intervall * 20000.) # this should be slightly less than nfev\n",
    "\n",
    "iXX1         = 2./5. * m1 * r1**2      # from the internet.\n",
    "iYY1         = iXX1\n",
    "iZZ1         = iXX1\n",
    "rhodtstreet1 = 1. # unimportant value, will be overwritten before used.\n",
    "\n",
    "pL_vals1  = [m1, mo1, 9.81, r1, amplitude1, frequenz1, alpha1, beta1, gamma1, iXX1, iYY1, iZZ1\n",
    "           ] + [nyb1, nys1, Eb1, Es1, ctau1, mu1] + [rhodtstreet1]\n",
    "pL1_vals1 = [m1, mo1, 9.81, r1, amplitude1, frequenz1, alpha1, beta1, gamma1, iXX1, iYY1, iZZ1\n",
    "           ] + [nyb1, nys1, Eb1, Es1, ctau1, mu1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ball must contact the street at *one* point only. Here I check, that the radius of the ball, $r_1$ is smaller than the smallest osculation circle (Schmiegekreis) of the street.\\\n",
    "if the radius of the ball is close to the osculating radius, things 'break down'. I suppose numerical issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if frequenz1 > 0.:\n",
    "    def func2(x, args):\n",
    "# just needed to get the arguments matching for minimuze\n",
    "        return np.abs(r_max_lam(*x, *args)[0])\n",
    "\n",
    "    def func3(x, args):\n",
    "# just needed to get the arguments matching for minimuze\n",
    "        return np.abs(r_max_lam(*x, *args)[1])\n",
    "\n",
    "    x0 = (0.1, 0.1)      # initial guess\n",
    "    minimal1 = minimize(func2, x0, pL1_vals1)\n",
    "    minimal2 = minimize(func3, x0, pL1_vals1)\n",
    "\n",
    "    minimal = min(minimal1.get('fun'), minimal2.get('fun'))\n",
    "\n",
    "    if pL_vals1[3] < minimal:\n",
    "        print('selected radius = {} is less than minimally admissible radius = {:.4f}, hence o.k.'\n",
    "          .format(pL_vals1[3], minimal), '\\n')\n",
    "    else:\n",
    "        print('selected radius {} is larger than admissible radius {:.4f}, hence NOT o.k.'\n",
    "          .format(pL_vals1[3], minimal), '\\n')\n",
    "        raise Exception('Radius of wheel is too large')\n",
    "\n",
    "else:\n",
    "    print('the street is a horizontal plane')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial values of $x_{st}, z_{st}$ are found numerically. I iterate twice, to get improved results, sometimes.\\\n",
    "If the ball is not above the street, an exception is raised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func1(x0, args):\n",
    "    return loesung_lam(*x0, *args).reshape(2)\n",
    "\n",
    "x0 = [1., 1.]                                 # initial guess\n",
    "args = [xb1, yb1, zb1, amplitude1, frequenz1]\n",
    "for _ in range(2):\n",
    "    antwort = fsolve(func1, x0, args)\n",
    "    x0 = antwort                              # improved guess\n",
    "\n",
    "antwort = [antwort[0], antwort[1]]\n",
    "\n",
    "min_hoehe = gesamt_lam(xb1, zb1 , amplitude1, frequenz1)\n",
    "if yb1 > min_hoehe:\n",
    "    print('minimal height of the ball above P0 must be {:.3f} < yb1 = {:.3f}, so this is fine.'\n",
    "      .format(min_hoehe, yb1))\n",
    "else:\n",
    "    raise Exception('yb1 too small, must be at least {:.3f}'.format(min_hoehe))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numerical Integration**\n",
    "\n",
    "*Comments*\n",
    "- max_step must be small in *solve_ivp*, else the integration may miss the points of the collision of the ball with the street.\n",
    "- The values of the collision coordinate, $x_{st}, z_{st}$ are available only during integration. They are collected in *zusatz1*.\n",
    "- I iterate over the numerical solution for $x_{st}, z_{st}$ a couple of times, this sometimes may improve the accuracy of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "start2 = time.time()\n",
    "max_step = 0.0001\n",
    "pL_vals  = antwort + pL_vals1\n",
    "pL1_vals = antwort + pL1_vals1\n",
    "print('starting parameters are', '\\n', ['{:.3f}'.format(pL_vals[i]) for i in range(len(pL_vals))], '\\n')\n",
    "\n",
    "times = np.linspace(0, intervall, schritte)\n",
    "t_span = (0., intervall)\n",
    "\n",
    "y0 = [xb1, yb1, zb1, q11, q21, q31] + [uxb1, uyb1, uzb1, u11, u21, u31]\n",
    "print('initial values are', '\\n', ['{:.3f}'.format(y0[i]) for i in range(len(y0))], '\\n')\n",
    "\n",
    "zusatz1 = []\n",
    "x0 = antwort\n",
    "\n",
    "def gradient(t, y, args):\n",
    "    global x0\n",
    "    args1 = [y[0], y[1], y[2], args[6], args[7]]\n",
    "# Iteration, to, hopefully, improve the accuracy of the solution\n",
    "    for _ in range(2):\n",
    "        antwort = fsolve(func1, x0, args1)\n",
    "        x0 = antwort\n",
    "    args[0] = antwort[0]\n",
    "    args[1] = antwort[1]\n",
    "    zusatz1.append([t, args[0], args[1]])\n",
    "    \n",
    "    if max_step < abstand1_lam(*y, *args) - r1 < 2. * max_step:\n",
    "        args[-1] = rhodtstreet1_lam(*y, *pL1_vals)\n",
    "    sol = np.linalg.solve(MM_lam(*y, *args), force_lam(*y, *args))\n",
    "    return np.array(sol).T[0]\n",
    "        \n",
    "resultat1 = solve_ivp(gradient, t_span, y0, t_eval = times, args=(pL_vals,), method='BDF', max_step=max_step\n",
    "            , atol=1.e-9, rtol=1.e-9)\n",
    "\n",
    "resultat = resultat1.y.T\n",
    "print('resultat shape', resultat.shape, '\\n')\n",
    "event_dict = {-1: 'Integration failed', 0: 'Integration finished successfully', 1: 'some termination event'}\n",
    "print(event_dict[resultat1.status], ', message is:', resultat1.message, '\\n')\n",
    "\n",
    "print(\"To numerically integrate an intervall of {} sec the routine cycled {} times and it took {:.3f} sec \"\n",
    "      .format(intervall, resultat1.nfev, time.time() - start2), '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot any **generalized coordinates or speeds** you want to see.\\\n",
    "*times* contains many points in time, needed to get the energies right, and also the hysteresis curves. Here this is not needed, so I reduce the points considered to around *zeitpunkte*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# reduce the number of points of time to around zeitpunkte\n",
    "times2 = []\n",
    "resultat2 = []\n",
    "index2 = []\n",
    "\n",
    "#=======================\n",
    "zeitpunkte = 500\n",
    "#=======================\n",
    "\n",
    "reduction = max(1, int(len(times)/zeitpunkte))\n",
    "\n",
    "for i in range(len(times)):\n",
    "    if i % reduction == 0:\n",
    "        times2.append(times[i])\n",
    "        resultat2. append(resultat[i])\n",
    "schritte2 = len(times2)\n",
    "resultat2 = np.array(resultat2)\n",
    "times2 = np.array(times2)\n",
    "\n",
    "\n",
    "bezeichnung = ['xb1', 'yb1', 'zb1', 'q11', 'q21', 'q31'] + ['uxb1', 'uyb1', 'uzb1', 'u11', 'u21', 'u31']\n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "for i in (0, 1, 2, 9, 10, 11):\n",
    "    ax.plot(times2, resultat2[:, i], label=bezeichnung[i])\n",
    "ax.set_xlabel('time (sec)')\n",
    "ax.set_title('Generalized coordinates and / or speeds ')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The locations of the contact point, $x_{st}$, $z_{st}$ are available only during the integration, where I collect them in *zusatz1*, see the function *gradient* above.\\\n",
    "Here I match them as closely as I can to the points in time stored in *times*. *index* will store the locations, where the match is close.\\\n",
    "This takes a may take a long time to run, as there are essentially $\\frac{1}{2} \\cdot |times|^2$ operations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zusatz1 = np.array(zusatz1)\n",
    "zaehler = 0\n",
    "index = []\n",
    "\n",
    "for zeit in times:\n",
    "    zaehler1 = np.min(np.argwhere(np.array(zusatz1[zaehler :, 0] >= zeit)))\n",
    "    index.append(zaehler+zaehler1)\n",
    "    zaehler = zaehler1\n",
    "if len(index) != len(times[: resultat.shape[0]]):\n",
    "    raise Exception('Something happened')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the **energies** of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pot_np    = np.empty(resultat.shape[0])\n",
    "kin_np    = np.empty(resultat.shape[0])\n",
    "spring_np = np.empty(resultat.shape[0])\n",
    "total_np  = np.empty(resultat.shape[0])\n",
    "\n",
    "for i in range(resultat.shape[0]):\n",
    "    pL_vals[0]   = zusatz1[index[i], 1]\n",
    "    pL_vals[1]   = zusatz1[index[i], 2]\n",
    "    pot_np[i]    = pot_lam(*[resultat[i, j] for j in range(12)], *pL_vals)\n",
    "    kin_np[i]    = kin_lam(*[resultat[i, j] for j in range(12)], *pL_vals)\n",
    "    spring_np[i] = spring_lam(*[resultat[i, j] for j in range(12)], *pL_vals)\n",
    "    total_np[i]  = pot_np[i] + kin_np[i] + spring_np[i]\n",
    "\n",
    "if ctau1 == 1. and mu1 == 0.:\n",
    "    max_total, min_total = np.max(total_np), np.min(total_np)\n",
    "    print('Max. deviation of the total energy from being constant is {:.2e} % of max. total energy'.\n",
    "          format((max_total - min_total)/max_total * 100))\n",
    "    \n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "ax.plot(times[: resultat.shape[0]], pot_np, label='potential energy')\n",
    "ax.plot(times[: resultat.shape[0]], kin_np, label='kinetic energy')\n",
    "ax.plot(times[: resultat.shape[0]], spring_np, label='spring energy')\n",
    "ax.plot(times[: resultat.shape[0]], total_np, label='total energy')\n",
    "ax.set_xlabel('time (sec)')\n",
    "ax.set_title('Energies of the system')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I plot the speed component of the ball in the $\\overline{Pst \\to Dmc}$ direction.\\\n",
    "Obviously it jumps almost instantaneoulsly during impact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "test_np = np.empty(resultat.shape[0])\n",
    "for i in range(resultat.shape[0]):\n",
    "    pL1_vals[0]   = zusatz1[index[i], 1]\n",
    "    pL1_vals[1]   = zusatz1[index[i], 2]\n",
    "    test_np[i] = rhodtstreet1_lam(*[resultat[i, j] for j in range(resultat.shape[1])], *pL1_vals)\n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "ax.plot(times[: resultat.shape[0]], test_np)\n",
    "ax.set_xlabel('time (sec)')\n",
    "ax.set_title('speed component in impact direction');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the **hysteresis curves** of successive impacts of the ball with the street.\\\n",
    "The red numbers are the times, approximately, when the impacts took place.\\\n",
    "As there may be many impacts, the picture may look messy. With the parameter *hyst* one may set the number of curves to be plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#==========================\n",
    "hyst = 2\n",
    "#==========================\n",
    "\n",
    "#determine R2 of the a.m. formula    \n",
    "r_xst   = -(sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(xst))**2\n",
    "               )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(xst, 2)\n",
    "r_zst   = -(sm.S(1.) + (gesamt(xst, zst, amplitude, frequenz).diff(zst))**2 \n",
    "               )**sm.S(3/2)/gesamt(xst, zst, amplitude, frequenz).diff(zst, 2)\n",
    "\n",
    "r_xst_lam = sm.lambdify([xst, zst, amplitude, frequenz], r_xst, cse=True)\n",
    "r_zst_lam = sm.lambdify([xst, zst, amplitude, frequenz], r_zst, cse=True)\n",
    "\n",
    "   \n",
    "sigma_b = (1. - nyb1**2) / Eb1\n",
    "sigma_s = (1. - nys1**2) / Es1\n",
    "\n",
    "HC_kraft = []\n",
    "HC_displ = []\n",
    "HC_times = []\n",
    "zaehler  = 0\n",
    "zaehler1 = 0\n",
    "i0 = 0\n",
    "    \n",
    "for i in range(resultat.shape[0]):\n",
    "    pL_vals[0]    = zusatz1[index[i], 1]\n",
    "    pL_vals[1]    = zusatz1[index[i], 2]\n",
    "    pL1_vals[0]   = zusatz1[index[i], 1]\n",
    "    pL1_vals[1]   = zusatz1[index[i], 2]\n",
    "    \n",
    "    xst1 = r_xst_lam(pL_vals[0], pL_vals[1], amplitude1, frequenz1)\n",
    "    zst1 = r_zst_lam(pL_vals[0], pL_vals[1], amplitude1, frequenz1)\n",
    "    R2   = 0.5 * (xst1 + zst1)\n",
    "    wurzel = sm.sqrt(sm.Abs(r1 * R2 /(r1 + R2))) # theoretically no need for sm.Abs(...)\n",
    "    k0      = 4. / (3.*(sigma_b + sigma_s)) * wurzel\n",
    "    \n",
    "    abstand = r1  - abstand1_lam(*[resultat[i, j] for j in range(resultat.shape[1])], \n",
    "            *pL_vals)\n",
    "    if abstand < 0.:\n",
    "        i0 = i+1\n",
    "\n",
    "    if abstand >= 0. and i0 == i:\n",
    "        zaehler1 += 1\n",
    "        walldt = rhodtstreet1_lam(*[resultat[i, j] for j in range(resultat.shape[1])], *pL1_vals)\n",
    "    try:\n",
    "        if zaehler1 > hyst:\n",
    "            raise Rausspringen\n",
    "        if abstand >= 0.:\n",
    "            rhodt = rhodtstreet1_lam(*[resultat[i, j] for j in range(resultat.shape[1])], *pL1_vals)\n",
    "            kraft0 = k0 * abstand**(3/2) * (1. + 3./2. * (1 - ctau1) * rhodt/walldt)\n",
    "            HC_displ.append(abstand)\n",
    "            HC_kraft.append(kraft0)\n",
    "            HC_times.append((zaehler, times[i]))\n",
    "            zaehler +=1\n",
    "    except Rausspringen:\n",
    "        break\n",
    "    \n",
    "HC_displ = np.array(HC_displ)\n",
    "HC_kraft = np.array(HC_kraft)\n",
    "\n",
    "# print only, if there were collisions with the street.    \n",
    "if len(HC_displ) != 0:\n",
    "    fig, ax = plt.subplots(figsize=(10,5))\n",
    "    ax.plot(HC_displ, HC_kraft, color='green')\n",
    "    ax.set_xlabel('penetration depth (m)')\n",
    "    ax.set_ylabel('contact force (N)')\n",
    "    ax.set_title(f'hysteresis curves of successive impacts of the ball_1 with the street, ctau = {ctau1}, mu = {mu1}')\n",
    "    \n",
    "    zeitpunkte = 8\n",
    "    reduction = max(1, int(len(HC_times)/zeitpunkte))\n",
    "    for k in range(len(HC_times)):\n",
    "        if k % reduction == 0:\n",
    "            coord  = HC_times[k][0]\n",
    "            ax.text(HC_displ[coord], HC_kraft[coord], f'{HC_times[k][1]:.4f}', color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Animate the movement*\\\n",
    "The color of the ball indicates its height above $P_0$\\\n",
    "*times* may have many entries, which makes the animation, at least if *HTML(..)* has to be used, take a long time to calculate. I reduce the number of entries to around *zeitpunkte*, maybe at the risk of loosing some accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "CP_pos_lam = sm.lambdify(qL + pL, [me.dot(Dmc.pos_from(P0), uv) for uv in N], cse=True)\n",
    "Po_pos_lam = sm.lambdify(qL + pL, [me.dot(Po.pos_from(P0), uv) for uv in (N.x, N.z)], cse=True)\n",
    "\n",
    "# reduce the number of points of time to around zeitpunkte\n",
    "times2 = []\n",
    "resultat2 = []\n",
    "index2 = []\n",
    "\n",
    "#=======================\n",
    "zeitpunkte = 500\n",
    "#=======================\n",
    "\n",
    "reduction = max(1, int(len(times)/zeitpunkte))\n",
    "\n",
    "for i in range(len(times)):\n",
    "    if i % reduction == 0:\n",
    "        times2.append(times[i])\n",
    "        resultat2. append(resultat[i])\n",
    "        index2.append(index[i]) \n",
    "schritte2 = len(times2)\n",
    "resultat2 = np.array(resultat2)\n",
    "times2 = np.array(times2)\n",
    "\n",
    "CPx  = np.empty(schritte2)\n",
    "CPy  = np.empty(schritte2)\n",
    "CPy1 = np.empty(schritte2)\n",
    "CPz  = np.empty(schritte2)\n",
    "\n",
    "Pox = np.empty(schritte2)\n",
    "Poz = np.empty(schritte2)\n",
    "\n",
    "for l in range(schritte2):\n",
    "    pL_vals[0]   = zusatz1[index2[l], 1]\n",
    "    pL_vals[1]   = zusatz1[index2[l], 2]\n",
    "    CPx[l] = CP_pos_lam(*[resultat2[l, j] for j in range(resultat2.shape[1])], *pL_vals)[0]\n",
    "    CPy[l] = CP_pos_lam(*[resultat2[l, j] for j in range(resultat2.shape[1])], *pL_vals)[1]\n",
    "    CPz[l] = CP_pos_lam(*[resultat2[l, j] for j in range(resultat2.shape[1])], *pL_vals)[2]\n",
    "    \n",
    "    Pox[l] = Po_pos_lam(*[resultat2[l, j] for j in range(resultat2.shape[1])], *pL_vals)[0]\n",
    "    Poz[l] = Po_pos_lam(*[resultat2[l, j] for j in range(resultat2.shape[1])], *pL_vals)[1]\n",
    "    \n",
    "\n",
    "# needed to give the picture the right size.\n",
    "xmin1 = min(CPx)\n",
    "xmin2 = min(CPz)\n",
    "xmin = min(xmin1, xmin2)\n",
    "\n",
    "xmax1 = max(CPx)\n",
    "xmax2 = max(CPz)\n",
    "xmax = max(xmax1, xmax2)\n",
    "\n",
    "for i in range(len(CPy)):\n",
    "    CPy1[i] = int(CPy[i] * 100.)\n",
    "ymin = min(CPy1)\n",
    "ymax = max(CPy1)\n",
    "\n",
    "# This is to asign colors of 'plasma' to the points.\n",
    "Test = mp.colors.Normalize(ymin, ymax)\n",
    "Farbe = mp.cm.ScalarMappable(Test, cmap='plasma')\n",
    "farbe1 = Farbe.to_rgba(CPy1[0]*100)    # color of the starting position\n",
    "    \n",
    "def animate_pendulum(zeit, x1, y1, z1, Pox1, Poz1):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={'aspect': 'equal'})\n",
    "    fig.colorbar(Farbe, label='Height of the ball above ground level \\n Factor = 100, the real height are the numbers divided by factor', shrink=0.9, ax=ax)\n",
    "    ax.axis('on')\n",
    "    ax.set_xlim(xmin - 2., xmax + 2.)\n",
    "    ax.set_ylim(xmin - 2., xmax + 2.)\n",
    "    ax.set_xlabel(' X Axis', fontsize=18)\n",
    "    ax.set_ylabel('Z Axis', fontsize=18)\n",
    "\n",
    "    line1, = ax.plot([], [], 'o', color=farbe1, markersize=35)                     # ball\n",
    "    line2, = ax.plot([], [], color='blue', linewidth=0.25)           # to trace the movement of CP\n",
    "    line3, = ax.plot([], [], 'o', color='red', markersize=5)          # the observer\n",
    "    \n",
    "    def animate(i):\n",
    "        farbe2 = Farbe.to_rgba(y1[i]*100)                                # color of the actual point at time i\n",
    "        ax.set_title('running time {:.2f} sec, with rumpel = {}'.format(zeit[i], rumpel), fontsize=15)\n",
    "        line1.set_data([x1[i]], [z1[i]])\n",
    "        line1.set_color(farbe2)\n",
    "        line2.set_data(x1[: i], z1[: i])  \n",
    "        line3.set_data([Pox1[i]], [Poz1[i]])\n",
    "        \n",
    "        return line1, line2, line3, \n",
    "\n",
    "    anim = animation.FuncAnimation(fig, animate, frames=len(times2),\n",
    "                                   interval=1000*max(times2) / len(times2),\n",
    "                                   blit=True)\n",
    "    plt.close(fig)\n",
    "    return anim\n",
    "\n",
    "anim = animate_pendulum(times2, CPx, CPy, CPz, Pox, Poz)\n",
    "print('It took {:.3f} sec, before HTML, to run the program'.format(time.time() - start))\n",
    "HTML(anim.to_jshtml())    # needed, when run on an iPad, I know no other way to do it. It is SLOW!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}