{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# billiard trajectory in a circle (animated)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nfrom collections import deque\n\ndef df(x):\n    return np.array([2*x[0], 2*x[1]])\n\ndef intersection(vec0, coord):\n    a = vec0[1]/vec0[0]\n    if np.abs(a) > 1e7:\n        z = np.array([x1, -y0])\n\n    x0 = coord[0]\n    y0 = coord[1]\n    sg = np.sign(vec0[0])\n    D = 1 - a**2*(-1+x0**2)+ 2*a*x0*y0 -y0**2\n    x1 = - ( -a**2*x0 + a*y0 - sg * np.sqrt(D) ) / (1 + a**2)\n\n    y1 = (y0 -a * ( x0 - sg* np.sqrt(D) ) )/ (1 + a**2)\n    z = np.array([x1, y1])\n    return z\n\nfig = plt.figure(figsize=(5, 5))\nax = fig.add_subplot(autoscale_on=False)\nax.set_aspect('equal')\nax.axis('off')\n\ntheta = np.linspace(-np.pi, np.pi, 100)\nz = np.exp(1j*theta)\nax.plot(z.real, z.imag, '-k',lw=3)\nax.set_xlim(-1.1, 1.1)\nax.set_ylim(-1.1, 1.1)\n\npoint0, = ax.plot([], [], 'o')\ntrace0, = ax.plot([], [], '-', lw=2, ms=2, c=point0.get_color())\n\n\ndt = 0.2\n\nhistory_len = 50\nhistory_x0, history_y0 = deque(maxlen=history_len), deque(maxlen=history_len)\n\nx0 = np.array([0.1, 0.4])\nvec0 = np.array([0.1, 0.02])\n\n\ntraj0 = [np.array([x0[0]]), np.array([x0[1]])]\n\nfor i in range(300):\n    x1 = x0 + vec0*dt\n\n    if x1[0]**2 + x1[1]**2 > 1:\n        x1 = intersection(vec0, x0)\n        n = -df(x1)\n        pn = -np.dot(vec0, n)/np.dot(n, n) * n\n        vec0 = vec0 + 2*pn\n\n    traj0[0] = np.append(traj0[0], x1[0])\n    traj0[1] = np.append(traj0[1], x1[1])\n    x0 = x1.copy()\n\n\ndef animate(i):\n    traj_x0,traj_y0 = traj0[0], traj0[1]\n\n    thisx0 = traj_x0[i]\n    thisy0 = traj_y0[i]\n\n\n    if i == 0:\n        history_x0.clear()\n        history_y0.clear()\n\n\n    history_x0.appendleft(thisx0)\n    history_y0.appendleft(thisy0)\n\n    point0.set_data(thisx0, thisy0)\n    trace0.set_data(history_x0, history_y0)\n\n    return point0, trace0\n\nani = animation.FuncAnimation(\n    fig, animate, len(traj0[0]), interval=200*dt, blit=True) #, len(y1), interval=dt*1000, blit=True)\n\nplt.show()"
      ]
    }
  ],
  "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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}