{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# parabolic reflector (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\n\n\nfig, ax = plt.subplots(1,1)\n\n\ndef y(x):\n    return x**2/2\n\ndef df(x,y):\n    return np.array([-x, 1])\n\n\nsample=10\nplot_l = []\nhist = []\nhist_len = 40\nax.axis('off')\nfig.patch.set_alpha(0.)\nax.plot(0, 0.5, 'ok', ms=10)\nfor i in range(sample):\n    point, = ax.plot([], [], '-', ms=10)\n    trace, = ax.plot([], [], '-', lw=2, c=point.get_color())\n    hist_x, hist_y = deque(maxlen=hist_len), deque(maxlen=hist_len)\n    plot_l.append(point)\n    plot_l.append(trace)\n    hist.append(hist_x)\n    hist.append(hist_y)\n\n\ndt = 0.02\n\nx0 = np.linspace(-1.5, 1.5, 10)\ny0 = np.array([2]*sample)\np = np.array([x0,y0])\n\ntraj_x = [np.array([])]*sample\ntraj_y = [np.array([])]*sample\n\nxx = np.linspace(-2, 2, 100)\nax.plot(xx, y(xx), '-k',lw=3)\nax.text(1.7,1.2, r\"$y=\\frac{x^2}{2}$\",fontsize=15)\n\nfor i, z in enumerate(zip(x0, y0)):\n    vec = np.array([0,-1])\n    for t in range(180):\n        z = z + dt*vec\n        x0, y0 = z\n        if y0 < y(x0):\n            x0 = z[0]\n            y0 = y(x0)\n            n = df(x0, y0)\n            pn = - np.dot(vec, n)/np.dot(n, n) * n\n            vec = vec + 2*pn\n        if  np.all(np.abs(np.array([x0,y0]) - np.array([0,0.5])) < 0.02):\n            vec = np.array([0,0])\n\n        traj_x[i] = np.append(traj_x[i], x0)\n        traj_y[i] = np.append(traj_y[i], y0)\n\n\n\ndef animate(i):\n    for j in range(sample):\n        thisx = traj_x[j][i]\n        thisy = traj_y[j][i]\n\n        if i == 0:\n            hist[2*j].clear()\n            hist[2*j+1].clear()\n\n        hist[2*j].append(thisx)\n        hist[2*j+1].append(thisy)\n\n        plot_l[2*j].set_data(thisx, thisy)\n        plot_l[2*j+1].set_data(hist[2*j], hist[2*j+1])\n\nani = animation.FuncAnimation(fig, animate, len(traj_x[0]),\n                              interval=dt*1000,\n#                              blit=True\n                              )\n\"\"\"\nani.save('anim_parabolic.mp4',\n#         writer='pillow'\n         writer=\"ffmpeg\",\n#         codec=\"png\",\n         savefig_kwargs={\"transparent\": True, 'facecolor': 'none'}\n         )\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.8.20"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}