{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# (classical) circular billard\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\n\nfig, ax = plt.subplots(1,1, figsize=(6,6))\n\ntheta = np.linspace(-np.pi, np.pi, 100)\nz = np.exp(1j*theta)\nax.plot(z.real, z.imag, '-')\nax.set_xlim(-1,1)\nax.set_ylim(-1,1)\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\n#x0 = np.array([0, 1/np.sqrt(2)])\nx0 = np.array([0.5, 0.1])\nvec0 = np.array([-0.2, 0.2])\n\nini_x0 = x0.copy()\nini_vec0 = vec0.copy()\n\nx1 = intersection(vec0, x0)\n\ntraj = [np.array([]), np.array([])]\n\ntraj[0] = np.append(traj[0], x0[0])\ntraj[1] = np.append(traj[1], x0[1])\n\ntraj[0] = np.append(traj[0], x1[0])\ntraj[1] = np.append(traj[1], x1[1])\n\nn = -df(x1)\npn = -np.dot(vec0, n)/np.dot(n, n) * n\nvec1 = vec0 + 2*pn \n\n#ax.quiver(x1[0], x1[1], n[0], n[1], angles='xy', scale_units='xy', scale=10)\n#ax.quiver(x1[0], x1[1], vec1[0], vec1[1], angles='xy', scale_units='xy', scale=1)\n\nfor n in range(30):\n    x0 = x1.copy()\n    vec0 = vec1.copy()\n\n    x1 = intersection(vec0, x0)\n    traj[0] = np.append(traj[0], x1[0])\n    traj[1] = np.append(traj[1], x1[1])    \n\n    n = -df(x1)\n    pn = -np.dot(vec0, n)/np.dot(n, n) * n\n    vec1 = vec0 + 2*pn \n\n    xx = np.linspace(x0[0], x1[0], 2)\n    yy = vec0[1]/vec0[0]*(xx - x0[0]) + x0[1]\n\n\nax.plot(traj[0], traj[1],'-', alpha=0.7)\n\nax.plot(ini_x0[0], ini_x0[1], 'ok')\nax.quiver(ini_x0[0], ini_x0[1], ini_vec0[0], ini_vec0[1], angles='xy', scale_units='xy', scale=0.5)\n\n\nax.grid()\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
}