{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# rectangular 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\nx = np.array([-3/2, 3/2, 3/2, -3/2, -3/2])\ny = np.array([1, 1, -1, -1, 1])\n\nax.plot(x, y, '-k', lw=3)\nax.set_ylim(-2,2)\nax.set_xlim(-2,2)\n\ndef intersection(vec0, coord):\n    a = vec0[1]/vec0[0]\n    x0 = coord[0]\n    y0 = coord[1]\n    sg = np.sign(vec0[1])\n    x1 = - ( - sg - a*x0 + y0 ) / a\n    y1 = sg\n\n    if np.abs(x1) >  3/2:\n        sg = np.sign(vec0[0])\n        x1 = sg*3/2\n        y1 = (sg*3*a - 2*a*x0 + 2*y0)/2\n\n    return np.array([x1, y1])\n\ndef df(x0):\n    x = x0[0]\n    y = x0[1]\n\n    if np.abs(np.abs(y) - 1) < 1e-10:\n        z = np.array([0, -1])\n    elif np.abs(np.abs(x) - 3/2) < 1e-10:\n        z = np.array([-1, 0])\n\n    return z\n    \n\ntraj = [np.array([]), np.array([])]\n\nx0 = np.array([0.1, 0.1])\nvec0 = np.array([-np.pi/10, 0.1])\n\ninit_x = x0.copy()\ninit_vec = vec0.copy()\n\ntraj[0] = np.append(traj[0], x0[0])\ntraj[1] = np.append(traj[1], x0[1])\n\nx1 = intersection(vec0, x0)\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=1)\n#ax.quiver(x1[0], x1[1], vec1[0], vec1[1], angles='xy', scale_units='xy', scale=1)\n\nfor i in range(10):\n    x0 = x1.copy()\n    vec0 = vec1.copy()\n    x1 = intersection(vec0, x0)\n\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    \n    #ax.quiver(x1[0], x1[1], vec1[0], vec1[1], angles='xy', scale_units='xy', scale=1)\n\nax.plot(traj[0], traj[1], '-')\nax.plot(init_x[0], init_x[1], 'ok')\nax.quiver(init_x[0], init_x[1], init_vec[0], init_vec[1], angles='xy', scale_units='xy', scale=1)\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
}