{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# (classical) stadium billiard\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy import optimize\n\n\nfig, ax = plt.subplots(1,1, figsize=(6,6))\n\ntheta1 = np.linspace(-np.pi/2, np.pi/2, 100)\nz1 = np.exp(1j*theta1) + 1\n\ntheta2 = np.linspace(np.pi/2, 3*np.pi/2, 100)\nz2 = np.exp(-1j*theta2) - 1\n\nx = np.array([-2, 2, 2, -2, -2])\ny = np.array([1, 1, -1, -1, 1])\n\nax.plot(z1.real, z1.imag, '-k', lw=3)\nax.plot(z2.real, z2.imag, '-k', lw=3)\nax.plot(x, y, '-k',lw=3)\n\nax.set_ylim(-2,2)\nax.set_xlim(-2,2)\n\ndef semi_cf(z, a, c, x0, y0):\n    x,y = z\n    return [(x - c)**2 + y**2 - 1,\n            y - a*(x - x0) - y0]\n\ndef circ_df(z, c):\n    x, y = z\n    return np.array([ 2 * ( -c + x ), 2 * y])\n\ndef rect_intersection(vec0, coord):\n    a = vec0[1]/vec0[0]\n    x0,y0 = coord\n    sg = np.sign(vec0[1])\n    x1 = - ( - sg - a*x0 + y0 ) / a\n    y1 = sg\n    if x1 > 1.0:\n        f = lambda x: semi_cf(x, a, 1, x0, y0)\n        res = optimize.root(f, [x1, y1])\n        x1, y1 = res.x\n    elif x1 < -1.0:\n        f = lambda x: semi_cf(x, a, -1, x0, y0)\n        res = optimize.root(f, [x1, y1])\n        x1, y1 = res.x\n    return np.array([x1, y1])\n\ndef df_rect(x0):\n    x,y = x0\n    if np.abs(np.abs(y) - 1) < 1e-10:\n        z = np.array([0, -1])\n    elif np.abs(np.abs(x) - 2) < 1e-10:\n        z = np.array([-1, 0])\n    return z\n\ndef df(x0):\n    x,y = x0\n    if np.abs(x) < 1:\n        return df_rect(x0)\n    elif x > 0:\n        return circ_df(x0, 1)\n    elif x < 0:\n        return circ_df(x0, -1)\n    else:\n        print(\"something wrong\")\n\n\ndef intersection(vec0, coord):\n    x, y = coord\n    #if np.abs(x) < 1:\n    return rect_intersection(vec0, coord)\n\ntraj = [np.array([]), np.array([])]\n\nx0 = np.array([0.1, 0.1])\nvec0 = np.array([0.05, 0.15])\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\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=5)\n#ax.quiver(x1[0], x1[1], vec1[0], vec1[1], angles='xy', scale_units='xy', scale=1/2)\n\nfor i in range(30):\n    x0 = x1.copy()\n    vec0 = vec1.copy()\n\n    x1 = intersection(vec0, x0)\n    ax.plot(x1[0], x1[1], 'ok')\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    #ax.quiver(x1[0], x1[1], n[0], n[1], angles='xy', scale_units='xy', scale=5)\n    #ax.quiver(x1[0], x1[1], vec1[0], vec1[1], angles='xy', scale_units='xy', scale=1/2)\n\nax.plot(traj[0], traj[1], '-',alpha=0.4)\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/2)\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
}