{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ring Resonator Dynamics\n",
    "The following analysis is base on the model from Yariv's PTL, \"Critical Coupling and Its Control in Optical Waveguide-Ring Resonator Systems\", in 2002. We mainly focus on the dynamic of how it reach the static state. We start from the light reach the coupler at the bus wavguide.\n",
    "\n",
    "P.S.: We split the time to discrete step of each round trip of light in the ring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "T = 0.99\n",
    "K = 1 - T\n",
    "Alpha = 0.99\n",
    "t = np.sqrt(T) + 0.0j\n",
    "k = 0.0 + 1j * np.sqrt(K)\n",
    "alpha = np.sqrt(Alpha)\n",
    "theta = 2.0 * np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a function to generate the time series of a1, a2, b1 & a2:  \n",
    "**a[n] => [a1(n-step), a2(n-step)]**  \n",
    "**b[n] => [b1(n-step), b2(n-step)]**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18264 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cd14a583c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def propagate_in_cycle(t, k, alpha, theta, num_cycle = 1000):\n",
    "    a1_0 = 1.0 + 0.0j\n",
    "    a2_0 = 0.0 + 0.0j\n",
    "    a = np.zeros((num_cycle, 2), dtype=complex)\n",
    "    b = np.zeros((num_cycle, 2), dtype=complex)\n",
    "\n",
    "    a[0] = [a1_0, a2_0]\n",
    "    \n",
    "    coupler = np.array([[t, k],\n",
    "                        [-k.conjugate(), t.conjugate()]])\n",
    "\n",
    "    for n in range(num_cycle - 1):\n",
    "        b[n] = coupler@a[n]\n",
    "        a[n+1] = np.array([a1_0, b[n][1] * alpha * np.exp(1.j * theta)])\n",
    "\n",
    "    b[num_cycle - 1] = coupler@a[num_cycle - 1]\n",
    "    \n",
    "    return a, b\n",
    "\n",
    "a, b = propagate_in_cycle(t=t, k=k, alpha=alpha, theta=theta, num_cycle = 1000)\n",
    "plt.figure()\n",
    "plt.subplot(121)\n",
    "plt.plot(abs(a[:,0])**2, label='a1')\n",
    "\n",
    "plt.plot(abs(b[:,0])**2, label='b1')\n",
    "plt.yscale('log')\n",
    "plt.title('a1 & b1 in Bus WG')\n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "plt.plot(abs(a[:,1])**2, label='a2')\n",
    "plt.plot(abs(b[:,1])**2, label='b2')\n",
    "plt.title('a2 & b2 in Resonator')\n",
    "plt.yscale('linear')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can invitigate the change of bus waveguide output intensity & the ring intensity gain from the coupler by each cycle (normalized by the input intensity a1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18216 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cd15ca4c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "bus_output = abs(b[:, 0]/a[:,0])**2\n",
    "gain_per_cycle = (abs(b[:, 1])**2 - abs(a[:, 1])**2)/abs(a[:,0])**2\n",
    "\n",
    "plt.plot(bus_output, label='Bus WG Output')\n",
    "plt.plot(gain_per_cycle, label='Energy gain of ring after coupler per cycle')\n",
    "plt.plot(bus_output + gain_per_cycle, label=\"Sum of both\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the bus waveguide output power:  \n",
    "**enegry loss from the ring via coupler + enegry from a1 does not coupled to the ring**  \n",
    "\n",
    "is same as the enegry gain of the ring after coupler:  \n",
    "**enegry gain from bus waveguide + energy of the ring does not coupled out to bus waveguide**  \n",
    "\n",
    "Therefore, when the distructive interference happened in bus output port (b1) of the coupler, the energy will transfer to the port of the ring (b2). It must be true since the coupler is lossless (coupler martrix is unitary, i.e., |t|^2 + |k|^2 = 1)  \n",
    "It is an interesting constructive & distructive interference porporty of wave 😉  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Play with other parameters\n",
    "Enjoy~ 😊"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 17660 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cd15c2b780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "T2 = 0.999\n",
    "K2 = 1 - T2\n",
    "Alpha2 = 0.999\n",
    "t2 = np.sqrt(T2) + 0.0j\n",
    "k2 = 0.0 + 1j * np.sqrt(K2)\n",
    "alpha2 = np.sqrt(Alpha2)\n",
    "theta2 = 2.0 * np.pi\n",
    "a2, b2 = propagate_in_cycle(t=t2, k=k2, alpha=alpha2, theta=theta2, num_cycle = 10000)\n",
    "plt.figure()\n",
    "plt.subplot(121)\n",
    "plt.plot(abs(a2[:,0])**2, label='a1')\n",
    "\n",
    "plt.plot(abs(b2[:,0])**2, label='b1')\n",
    "plt.yscale('log')\n",
    "plt.title('a1 & b1 in Bus WG')\n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "plt.plot(abs(a2[:,1])**2, label='a2')\n",
    "plt.plot(abs(b2[:,1])**2, label='b2')\n",
    "plt.title('a2 & b2 in Resonator')\n",
    "plt.yscale('linear')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18428 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cd15b431d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "bus_output = abs(b2[:, 0]/a2[:,0])**2\n",
    "gain_per_cycle = (abs(b2[:, 1])**2 - abs(a2[:, 1])**2)/abs(a2[:,0])**2\n",
    "\n",
    "plt.plot(bus_output, label='Bus WG Output')\n",
    "plt.plot(gain_per_cycle, label='Energy gain of ring after coupler per cycle')\n",
    "plt.plot(bus_output + gain_per_cycle, label=\"Sum of both\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
