newtons-pendulum-Hunt-Crossley.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy.physics.mechanics as me\n",
    "import sympy as sm\n",
    "from scipy.integrate import odeint, solve_ivp\n",
    "from scipy.interpolate import pchip\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML\n",
    "import matplotlib\n",
    "import time\n",
    "matplotlib.rcParams['animation.embed_limit'] = 2**128"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is to model **Newton's Pendulum**. If the leftmost pedulum\n",
    "is moved out and released, it hits the pendulums resting, and only the rightmost one should swing, and\n",
    "back again.\n",
    "\n",
    "The pendulums which are hit do not remain totally motionless as the link below shows.\n",
    "\n",
    "https://en.wikipedia.org/wiki/File:Newton%E2%80%99s_cradle_slo_mo.webm\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 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",
    "\n",
    "As per the article, $n = \\frac{3}{3}$ 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",
    "**Notes**  \n",
    "1.\\\n",
    "$c_\\tau = 1.$ gives **Hertz's** solution to the impact problem, also described in the article\\\n",
    "2.\\\n",
    "For $n > 2$ only $c_\\tau = 1.$ gives satisfactory results. I do not know, why.\n",
    "\n",
    "**Variables**\n",
    "\n",
    "- $n$: number of pendulums\n",
    "- $q_1...q_n$: generalized coordinate for the pendulums\n",
    "- $u_1...u_n$: the angular speeds\n",
    "- $N$: frame of inertia\n",
    "- $A_i$: body fixed frame of the i-th pendulum\n",
    "\n",
    "- $m_i$: mass of the i-th pendulum\n",
    "- $r_0, l_0, k_0$: radius of pendulum, distance of the center of the pendulum to its suspension point, that is | $\\overline{P_{unten_i}P_{oben_i}}$ |, modulus of elasticity of the pendulum body\n",
    "- $i_{ZZ_i}$: moment of ineratia of the i-th pendulum\n",
    "- $P_{oben_i}$: suspension point of i-th pendulum\n",
    "- $P_{unten_i}$: center of mass of i-th pendulum\n",
    "\n",
    "\n",
    "- $c_\\tau$: the experimental constant needed for Hunt-Crossley\n",
    "- $rhodt_{max}$: the collission speed, to be determined during integration, needed for Hunt_Crossley "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HC_force(N, P1, P2, r, ctau, rhodtmax, k0):\n",
    "    '''\n",
    "Calculates the contact force exerted by P1 on P2, according to the Hunt-Crossley formula given above.\n",
    "I assume, that the contact force always acts along the line P1/P2. I think, this is a fair assymption if\n",
    "the colliding balls are homogenious.\n",
    "\n",
    "The variables in the list are\n",
    "\n",
    "- N is an me.ReferenceFrame, the inertial frame\n",
    "- P1, P2 are me.Point objects. They are the centers of two me.RigidBody objects, here assumed to be two\n",
    "  ball each of radius r\n",
    "- radius of the ball\n",
    "- ctau, the experimental constant needed\n",
    "- rhodtmax, the relative speeds of P1 to P2, right at the impact time, measured in N. \n",
    "  This has to be calculated numerically during the integration, I found no ther way.\n",
    "- k0 is the force constant\n",
    "\n",
    "When I take the magnitude of the penetration speed, called geschw, I 'loose' its direction. \n",
    "I use vorzeichen to get it.\n",
    "\n",
    "    '''\n",
    "    vektorP1P2 = P2.pos_from(P1)\n",
    "    abstand = vektorP1P2.magnitude() \n",
    "    richtung = vektorP1P2.normalize()        \n",
    "    \n",
    "    rho = 2. * r - abstand       # penetration. Positive if the two balls are in collision\n",
    "    geschw = vektorP1P2.diff(t, N)\n",
    "    vorzeichen = sm.sign(me.dot(geschw, vektorP1P2))\n",
    "    rhodt = vorzeichen * geschw.magnitude()\n",
    "    rho = sm.Max(rho, sm.S(0))   # if rho < 0., rho**(3/2) will give problems\n",
    "    \n",
    "    kraft = k0 * rho**(3/2) * (1. + 3./2. * (1 - ctau) * rhodt/rhodtmax\n",
    "        ) * sm.Heaviside(2. * r - abstand, 0.) * richtung\n",
    "    \n",
    "    return kraft   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "start = time.time()\n",
    "#-----------------------------------------------------------------------------------------------------------\n",
    "n = 2                      # number of pendulums, must the >=2\n",
    "#------------------------------------------------------------------------------------------------------------------------------------\n",
    "N = me.ReferenceFrame('N')\n",
    "P0 = me.Point('P0')\n",
    "P0.set_vel(N, 0)\n",
    "\n",
    "q = []\n",
    "u = []\n",
    "\n",
    "masse = []                                      # different particles may have different masses.\n",
    "inertia = []\n",
    "A = []\n",
    "for i in range(n):                              # Note: dates for point i are on place i-1\n",
    "    q.append(me.dynamicsymbols('q'+str(i), real=True))     # Rotation of frame i relative to N\n",
    "    u.append(me.dynamicsymbols('u'+str(i), real=True))     # angular velocity \n",
    "    masse.append(sm.symbols('m' + str(i)))\n",
    "    inertia.append(sm.symbols('iZZ' + str(i)))\n",
    "    A.append(me.ReferenceFrame('A' + str(i)))\n",
    "    \n",
    "g, r0, l0, k0, ctau = sm.symbols('g, r0, l0, k0, ctau', real=True)    \n",
    "t = sm.symbols('t')\n",
    "\n",
    "P_oben = []                # attachment point of pendulum on top\n",
    "P_unten = []               # pendulum point\n",
    "BODY = []\n",
    "traegheit = []             # this holds the moments of inertia of the pemdulum bodies\n",
    "rhodtmax = []              # relative speeds P_i / P_i+1 right at the impact\n",
    "\n",
    "for i in range(n):\n",
    "    A[i].orient_axis(N, q[i], N.z)\n",
    "    A[i].set_ang_vel(N, u[i]*N.z)\n",
    "    P_oben.append(P0.locatenew('P0' + str(i), 2*i*r0*N.x))\n",
    "    P_unten.append(P_oben[i].locatenew('P1' + str(i), l0*A[i].y))\n",
    "    traegheit.append(me.inertia(A[i], 0., 0., inertia[i]))\n",
    "    rhodtmax.append(sm.symbols('rhodtmax' + str(i)))\n",
    "    \n",
    "    BODY.append( (me.RigidBody('P1'+str(i)+'a', P_unten[i], A[i], masse[i], (traegheit[i], P_unten[i]))) )\n",
    "\n",
    "#Kinemat. equations\n",
    "kd = [*[u[i] - q[i].diff(t) for i in range(n)]]\n",
    "\n",
    "# Impact forces. Leftmost and rightmost pendulae are different from the ones in the middle\n",
    "FL = [(P_unten[0], HC_force(N, P_unten[1], P_unten[0], r0, ctau, rhodtmax[0], k0))]\n",
    "for i in range(1, n-1):\n",
    "    print('da')\n",
    "    FL.append((P_unten[i], HC_force(N, P_unten[i-1], P_unten[i], r0, ctau, rhodtmax[i-1], k0) + \n",
    "               HC_force(N, P_unten[i+1], P_unten[i], r0, ctau, rhodtmax[i], k0)  ))\n",
    "FL.append((P_unten[n-1], HC_force(N, P_unten[n-2], P_unten[n-1], r0, ctau, rhodtmax[n-2], k0)))\n",
    "\n",
    "# gravitational forces\n",
    "for i in range(n):\n",
    "    FL.append((P_unten[i], -masse[i] * g * N.y))\n",
    "\n",
    "# Get the speed of P_i / P_i+1- needed for the HC method.\n",
    "rhodt_impact = []\n",
    "distanz1 = []\n",
    "for i in range(n-1):\n",
    "    abstand = P_unten[i+1].pos_from(P_unten[i])\n",
    "    distanz1.append(abstand.magnitude())\n",
    "    geschw = (abstand.diff(t, N)).subs({sm.Derivative(l, t): j for l, j in zip(q, u)})\n",
    "    vorzeichen = sm.sign(me.dot(abstand, geschw))\n",
    "    rhodt_impact.append(vorzeichen * geschw.magnitude())\n",
    "\n",
    "#Kane's equations\n",
    "KM = me.KanesMethod(N, q_ind=q, u_ind=u, kd_eqs=kd)\n",
    "(fr, frstar) = KM.kanes_equations(BODY, FL)\n",
    "\n",
    "MM = KM.mass_matrix_full\n",
    "print('MM DS', me.find_dynamicsymbols(MM))\n",
    "print('MM free symbols', MM.free_symbols)\n",
    "print('MM contains {} operations'.format(sum([MM[i, j].count_ops(visual=False) for i in range(MM.shape[0])\n",
    "        for j in range(MM.shape[1])])), '\\n')\n",
    "\n",
    "force = KM.forcing_full\n",
    "print('force DS', me.find_dynamicsymbols(force))\n",
    "print('force free symbols', force.free_symbols)\n",
    "print('force contains {} operations'.format(sum([force[i].count_ops(visual=False) \n",
    "        for i in range(force.shape[0]) ])), '\\n')\n",
    "\n",
    "P_unten_list = [[me.dot(P_unten[i].pos_from(P0), uv) for uv in (N.x, N.y)] for i in range(n)]\n",
    "P_oben_list = [[me.dot(P_oben[i].pos_from(P0), uv) for uv in (N.x, N.y)] for i in range(n)]\n",
    "\n",
    "# Always good to look at the energy of the system\n",
    "pot_energie = sum([masse[i] * g * me.dot(P_unten[i].pos_from(P0), N.y) for i in range(n)])\n",
    "kin_energie = sum([BODY[i].kinetic_energy(N) for i in range(n)])\n",
    "\n",
    "\n",
    "# I must use the sm.Max(, 0) method, as a**(5/2) gives problems when a < 0.\n",
    "spring_energie = 0.\n",
    "for i in range(n-1):\n",
    "    deltas = sm.Max(2.*r0 - distanz1[i], sm.S(0.))\n",
    "    deltas = deltas**(5/2)\n",
    "    spring_energie += k0 * 2./5. * deltas * sm.Heaviside(2.*r0 - distanz1[i], 0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Force during penetration*\\\n",
    "This is needed only for the plotting of the *hysteresis curves* below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "vektorP1P2 = P_unten[1].pos_from(P_unten[0])\n",
    "abstand = vektorP1P2.magnitude() \n",
    "richtung = vektorP1P2.normalize()        \n",
    "\n",
    "geschw = vektorP1P2.diff(t, N)\n",
    "vorzeichen = sm.sign(me.dot(vektorP1P2, geschw))\n",
    "rho = 2. * r0 - abstand       # penetration. Positive if the two balls are in collision\n",
    "rhodt = vorzeichen * geschw.magnitude()\n",
    "rho = sm.Max(rho, sm.S(0))   # if rho < 0., rho**(3/2) will give problems\n",
    "    \n",
    "kraft = (k0 * rho**(3/2) * (1. + 3./2. * (1 - ctau) * rhodt/rhodtmax[0]\n",
    "        ) * sm.Heaviside(2. * r0 - abstand, 0.)).subs({sm.Derivative(q[0], t): u[0],\n",
    "                sm.Derivative(q[1], t): u[1]})\n",
    "print('kraft dynamic symbols', me.find_dynamicsymbols(kraft, reference_frame=N))\n",
    "print('kraft free symbols   ', kraft.free_symbols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the *sympy functions* are converted to *numpy functions* so numerical calculations may be done. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Lambdification\n",
    "\n",
    "#def myHeaviside(x):             \n",
    "#needed so lambdify knows what Heaviside is. (Trick from Jason Moore) Not\n",
    "# needed anymore as Heaviside has been implemented into sympy.\n",
    "# I do not delete it, so I know how to do it if some other function is not implemented in sympy.\n",
    "#    return np.heaviside(x, 0.5)\n",
    "\n",
    "qL = q + u\n",
    "pL = masse + inertia + [g, r0, l0, k0, ctau] + rhodtmax\n",
    "pLL = masse + inertia + [g, r0, l0, k0, ctau]\n",
    "\n",
    "MM_lam = sm.lambdify(qL + pL, MM, cse=True) #,  [{'Heaviside': myHeaviside}, 'numpy'] )  # Das 'numpy' ist notwendig.\n",
    "force_lam = sm.lambdify(qL + pL, force, cse=True) #, [{'Heaviside': myHeaviside}, 'numpy'] )\n",
    "\n",
    "P_oben_lam = sm.lambdify(qL + pL, P_oben_list, cse=True)\n",
    "P_unten_lam = sm.lambdify(qL + pL, P_unten_list, 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",
    "rhodt_impact_lam = sm.lambdify(qL + pLL, rhodt_impact, cse=True)\n",
    "distanz1_lam = sm.lambdify(qL + pLL, distanz1, cse=True)\n",
    "\n",
    "kraft_lam = sm.lambdify(qL + pL, kraft, cse=True)\n",
    "\n",
    "print('it took {:.3f} sec sec to establish Kanes equations'.format(time.time() - start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numerical Integration**\n",
    "\n",
    "All is standard.\\\n",
    "In *solve_ivp* using *method = 'BDF'* seems to be better here than using *method = 'Radau'*.\\\n",
    "*rhodtmax* is found numerically during the integration. I believe, this is the most 'tricky' part of the whole program.\n",
    "\n",
    "The number of result points given by solve_ivp, this value is given by *schritte* must be very large to get a nice plot of the *hysteresis curves*\\\n",
    "If I use the Young modulus for steel, around $E_Y = 2. \\cdot 10^{11} \\frac{N}{m^2}$, I would mave to make max_steps still much smaller, resulting in very long integration times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Numerical integration\n",
    "start1 = time.time()\n",
    "\n",
    "# Input variables\n",
    "#-------------------------------------------------------------------\n",
    "r01 = 1.0           # radius of the pendulum\n",
    "l01 = 10.           # distance of particles from the ceiling\n",
    "m01 = 1.            # mass of a pendulum\n",
    "ctau1 = 0.8         # given in the article\n",
    "\n",
    "intervall = 20                          \n",
    "\n",
    "q01 = [np.pi for i in range(n)] # all particles hang down verticvally, except the first one\n",
    "u01 = [0. for i in range(n)]    # all have angular velocity = 0. \n",
    "q01[0] = 0.5 * np.pi            # first particle moved\n",
    "#-------------------------------------------------------------------\n",
    "if n > 2:\n",
    "    ctau1 = 1.       # for n > 2, only this case works well\n",
    "    \n",
    "max_step = 0.001\n",
    "\n",
    "rhodtmax1 = [1. + k for k in range(len(rhodtmax))] # of no importance, could be anything\n",
    "\n",
    "nu = 0.28                       # Poisson's ratio, from the internet\n",
    "EY = 2.e7                       # units: N/m^2, Young's modulus from the internet, around 2e11 for steel\n",
    "sigma = (1 - nu**2) / EY\n",
    "k01 = 4. / (3.* (sigma + sigma)) * np.sqrt(r01**2/(r01 + r01))\n",
    "print(f'k01 = {k01:.3e}')\n",
    "\n",
    "iZZa = 0.25 * m01 * r01**2      # from the internet\n",
    "\n",
    "schritte = 300000     \n",
    "times = np.linspace(0, intervall, schritte)\n",
    "y0 = q01 + u01\n",
    "\n",
    "pL_vals  = [m01 for i in range(n)] + [iZZa for i in range(n)] + [9.8, r01, l01, k01, ctau1] + rhodtmax1\n",
    "\n",
    "def gradient(t, y, args):\n",
    "    \n",
    "# here I find rhodtmax, the speed just before the impact\n",
    "    for i in range(n-1):\n",
    "        hilfss = distanz1_lam(*y, *[args[k] for k in range(2*n+5)])[i] - 2.*args[2*n+1]\n",
    "        if (0. < hilfss < 0.001 and y[n+1+i] != 0.) or (0 < hilfss < 0.001 and y[n+i] != 0.):\n",
    "            args[2 * n + 5 + i] = rhodt_impact_lam(*y, *[args[k] for k in range(int(2*n+5))])[i]\n",
    "                    \n",
    "    vals = np.concatenate((y, args))\n",
    "    sol = np.linalg.solve(MM_lam(*vals), force_lam(*vals))\n",
    "    return np.array(sol).T[0]\n",
    "\n",
    "    \n",
    "resultat1 = solve_ivp(gradient, (0., float(intervall)), y0, args=(pL_vals,), t_eval=times, method='BDF', \n",
    "            max_step=max_step, atol=1.e-8, rtol=1.e-8)\n",
    "resultat = resultat1.y.T\n",
    "event_dict = {-1: 'Integration failed', 0: 'Integration finished successfully', 1: 'some termination event'}\n",
    "print(event_dict[resultat1.status])\n",
    "\n",
    "print('to calculate an intervall of {:.0f} sec it took {} loops and {:.3f} sec running time'.\n",
    "      format(intervall, resultat1.nfev, time.time() - start1))\n",
    "#=================================================================================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Plot the coordinates*\\\n",
    "As the number of data points is very large, much larger than needed here, I reduce them to $N_2$ to speed up plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "N2 = 5000\n",
    "N1 = int(resultat.shape[0] / N2)\n",
    "times1 = []\n",
    "resultat1 = []\n",
    "for i in range(resultat.shape[0]):\n",
    "    if i % N1 == 0:\n",
    "        times1.append(times[i])\n",
    "        resultat1.append(resultat[i])\n",
    "resultat1 = np.array(resultat1)\n",
    "times1 = np.array(times1)\n",
    "\n",
    "# print generalized coordinates\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "for i in range(2*n):\n",
    "    ax.plot(times1[: resultat1.shape[0]], resultat1[:, i], label='Koordinate {}'.format(i))\n",
    "ax.set_title(f'Generalized coordinates and angular velocities, ctau = {ctau1}', fontsize = 12)\n",
    "ax.legend();\n",
    "print(pL_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Plot the energies of the system*\n",
    "\n",
    "In order to get fairly constant total energy, the values of $a_{tol}, r_{tol}$ must be set ( in solve_ivp ) smaller than standard.\\\n",
    "As the impact times are very short, they do not always show in the plot: they are not returned by *solve_ivp*.\\\n",
    "for $EY \\approx 10^7$ the impacts are visible.\n",
    "\n",
    "As the number of results is very large, much larger than needed here, I reduce it to $N_2$ data points. This speeds up the plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "N2 = 20000\n",
    "N1 = int(resultat.shape[0] / N2)\n",
    "times1 = []\n",
    "resultat1 = []\n",
    "for i in range(resultat.shape[0]):\n",
    "    if i % N1 == 0:\n",
    "        times1.append(times[i])\n",
    "        resultat1.append(resultat[i])\n",
    "resultat1 = np.array(resultat1)\n",
    "times1 = np.array(times1)\n",
    "\n",
    "# print energies\n",
    "pot_np = np.empty(resultat1.shape[0])\n",
    "kin_np = np.empty(resultat1.shape[0])\n",
    "spring_np = np.empty(resultat1.shape[0])\n",
    "total_np = np.empty(resultat1.shape[0])\n",
    "\n",
    "\n",
    "for i in range(resultat1.shape[0]):\n",
    "    pot_np[i] = pot_lam(*[resultat1[i, j]  for j in range(resultat1.shape[1])], *pL_vals)\n",
    "    kin_np[i] = kin_lam(*[resultat1[i, j]  for j in range(resultat1.shape[1])], *pL_vals)\n",
    "    spring_np[i] = spring_lam(*[resultat1[i, j]  for j in range(resultat1.shape[1])], *pL_vals)\n",
    "    total_np[i] = pot_np[i] + kin_np[i] + spring_np[i]\n",
    "\n",
    "min_total = np.abs(np.min(total_np))\n",
    "max_total = np.abs(np.max(total_np))\n",
    "\n",
    "print('Max. deviation of 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(times1[: resultat1.shape[0]], pot_np, label='potential energy')\n",
    "ax.plot(times1[: resultat1.shape[0]], kin_np, label='kinetic energy')\n",
    "ax.plot(times1[: resultat1.shape[0]], spring_np, label='spring energy')\n",
    "ax.plot(times1[: resultat1.shape[0]], total_np, label='total energy')\n",
    "ax.set_title('Energies of the system of {} bodies, ctau = {}'.format(n, ctau1), fontsize=12)\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hysteresis curve**\\\n",
    "Here the hysteresis curves of the first ball are plotted.\\\n",
    "They get 'smaller' at subsequet impacts, as is to be expected: Energy is lost at each impact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "HC_kraft = []\n",
    "HC_displ = []\n",
    "\n",
    "for i in range(resultat.shape[0]):\n",
    "    if distanz1_lam(*[resultat[i, j] for j in range(resultat.shape[1])], *[pL_vals[k]\n",
    "            for k in range(2*n+5)])[0] - 2.*pL_vals[2*n+1] <= 0.:\n",
    "        HC_displ.append(-(distanz1_lam(*[resultat[i, j] for j in range(resultat.shape[1])], \n",
    "            *[pL_vals[k] for k in range(2*n+5)])[0] - 2.*pL_vals[2*n+1]))\n",
    "        HC_kraft.append(kraft_lam(*[resultat[i, j] for j in range(resultat.shape[1])], *pL_vals))\n",
    "\n",
    "HC_displ = np.array(HC_displ)\n",
    "HC_kraft = np.array(HC_kraft)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "ax.plot(HC_displ, HC_kraft)\n",
    "ax.set_xlabel('penetration depth (m)')\n",
    "ax.set_ylabel('contact force (Nm)')\n",
    "ax.set_title('hysteresis curves of successive impacts, ctau = {}'.format(ctau1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Animation*\n",
    "\n",
    "I limit the number of time points to be considered to around *zeitpunkte*, do reduce the time it takes to plot. HTML is VERY slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Animation\n",
    "\n",
    "times2 = []\n",
    "resultat2 = []\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",
    "\n",
    "schritte2 = len(times2)\n",
    "resultat2 = np.array(resultat2)\n",
    "times2 = np.array(times2)\n",
    "print('number of points considered:',len(times2))\n",
    "\n",
    "Koordinaten_unten = [np.zeros((schritte2, 2))for i in range(n)]\n",
    "Koordinaten_oben = [np.zeros(2) for i in range(n)]\n",
    "\n",
    "for i in range(n):\n",
    "    for j in range(schritte2):\n",
    "        Koordinaten_unten[i][j, 0] = P_unten_lam(*[resultat2[j, k]     # X coordinate of i_th pendulum\n",
    "                for k in range(resultat.shape[1])], *pL_vals)[i][0]       \n",
    "        Koordinaten_unten[i][j, 1] = P_unten_lam(*[resultat2[j, k]     # Y coordinate of i_th pendulum\n",
    "                for k in range(resultat.shape[1])], *pL_vals)[i][1]\n",
    "        \n",
    "    Koordinaten_oben[i][0] = P_oben_lam(*[resultat2[0, k]   \n",
    "                for k in range(resultat.shape[1])], *pL_vals)[i][0]       \n",
    "    Koordinaten_oben[i][1] = P_oben_lam(*[resultat2[0, k]   \n",
    "                for k in range(resultat.shape[1])], *pL_vals)[i][1] \n",
    "\n",
    "obenx = np.zeros(n)\n",
    "obeny = np.zeros(n)\n",
    "for i in range(n):\n",
    "    obenx[i] = Koordinaten_oben[i][0]\n",
    "    obeny[i] = Koordinaten_oben[i][1]\n",
    "\n",
    "def animate_pendulum(times2, Koordinaten_unten):\n",
    "    xmin = min([Koordinaten_unten[0][j, 0] for j in range(schritte2)])\n",
    "    xmax = max([Koordinaten_unten[n-1][j, 0] for j in range(schritte2)])\n",
    "    \n",
    "    ymin = min([Koordinaten_unten[0][j, 1] for j in range(schritte2)])\n",
    "    ymax = max([Koordinaten_unten[n-1][j, 1] for j in range(schritte2)])\n",
    "    \n",
    "    kleinst = min(xmin, ymin)\n",
    "    groesst = max(xmax, ymax)\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(15, 15))\n",
    "    ax.set_aspect('equal')\n",
    "    ax.axis('on')\n",
    "    ax.set(xlim=(kleinst - 1., groesst + 1.), ylim=(kleinst - 1., 1.))\n",
    "    ax.scatter(obenx, obeny, marker='o', linewidths=5, color='black') # points on the ceiling\n",
    "    \n",
    "    LINE1 = []\n",
    "    LINE2 = []\n",
    "    for i in range(n):\n",
    "# picking the 'right' radius of the discs I do by trial and error. I did not try to get a formula\n",
    "        line1, = ax.plot([], [], 'o', markersize= 45 + 5.*(8-n))\n",
    "        line2, = ax.plot([], [], lw=0.5, color='blue')\n",
    "        LINE1.append(line1)\n",
    "        LINE2.append(line2)\n",
    "\n",
    "    def animate(i):\n",
    "        ax.set_title('Pendulum of {} bodies, running time {:.2f} sec'.format(n, i/schritte2 * intervall)\n",
    "                     , fontsize=20)\n",
    "        for j in range(n):\n",
    "            x = Koordinaten_unten[j][:, 0]\n",
    "            y = Koordinaten_unten[j][:, 1]\n",
    "            LINE1[j].set_data(x[i], y[i])\n",
    "            LINE2[j].set_data((obenx[j], x[i]), (obeny[j], y[i]))\n",
    "        return LINE1 + LINE2\n",
    "\n",
    "    anim = animation.FuncAnimation(fig, animate, frames=schritte2,\n",
    "                                   interval=1000*times2.max() / schritte2,\n",
    "                                   blit=True)\n",
    "    plt.close(fig)\n",
    "    return anim\n",
    "\n",
    "anim = animate_pendulum(times2, Koordinaten_unten)\n",
    "print(f'it took {time.time() - start :.3f} sec to run the program, BEFORE HTML')\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}