fiber

environment.yml

dependencies:
  - python=3.7
  - ipywidgets==7.5.1
  - bqplot==0.11.0
  - numpy== 1.19.2
  - scipy==1.5.2

fiber_modes_v01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fiber Modes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assign U to any electric or magnetic field, the equation to solved is called the Helmholtz equation:\n",
    "\n",
    "$$\n",
    "\\nabla^2 U + n^2(r)k_0^2 U=0\n",
    "$$\n",
    "\n",
    "These are the solutions:\n",
    "\n",
    "$U(r,\\phi,z)=u(r)e^{-jl\\phi}e^{-j\\beta z}$  for  $l=0,\\pm 1,\\pm 2...$\n",
    "\n",
    "$$\n",
    "u(r) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "        J_l(k_T r) & r<a  \\\\\n",
    "        K_l(\\gamma r) & r>a\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "where J and K are the Bessel functions of the first and second kind respectively.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-11T23:40:33.739957Z",
     "start_time": "2021-02-11T23:40:32.566952Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.special as spe\n",
    "from scipy import optimize\n",
    "\n",
    "import bqplot\n",
    "from bqplot import pyplot as bqplt\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameters and Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-11T23:40:33.771957Z",
     "start_time": "2021-02-11T23:40:33.742953Z"
    }
   },
   "outputs": [],
   "source": [
    "pi = np.pi\n",
    "\n",
    "xx = np.linspace(-1.7, 1.7, 60)\n",
    "yy = np.linspace(-1.7, 1.7, 60)\n",
    "\n",
    "x_mesh, y_mesh = np.meshgrid(xx, yy)\n",
    "r_mesh = np.sqrt(x_mesh ** 2 + y_mesh ** 2)\n",
    "phi_mesh = np.arctan2(y_mesh, x_mesh)\n",
    "\n",
    "ones_mesh = np.ones((len(xx), len(yy)))\n",
    "zeros_mesh = np.zeros((len(xx), len(yy)))\n",
    "in_core_mesh = ones_mesh.copy()\n",
    "in_core_mesh[r_mesh > 1] = zeros_mesh[r_mesh > 1]  # mask core with ones\n",
    "in_clad_mesh = ones_mesh.copy()\n",
    "in_clad_mesh[r_mesh <= 1] = zeros_mesh[r_mesh <= 1]  # mask cladding with ones\n",
    "\n",
    "# this is to plot the core prerimeter later\n",
    "phi_core_shape = np.linspace(0, 2 * pi, 60)\n",
    "x_core_shape = 1 * np.cos(phi_core_shape)\n",
    "y_core_shape = 1 * np.sin(phi_core_shape)\n",
    "\n",
    "n_cladding = 1.444"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Auxiliary Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-11T23:40:33.819957Z",
     "start_time": "2021-02-11T23:40:33.774957Z"
    }
   },
   "outputs": [],
   "source": [
    "def zero_func(X, V, L):\n",
    "    Y = np.sqrt(V ** 2 - X ** 2)\n",
    "    return X * spe.jv(L + 1, X) / spe.jv(L, X) - Y * spe.kv(L + 1, Y) / spe.kv(L, Y)\n",
    "\n",
    "\n",
    "def find_zeros_exact(X, Y, V, L):\n",
    "    f = X * spe.jv(L + 1, X) / spe.jv(L, X) - Y * spe.kv(L + 1, Y) / spe.kv(L, Y)\n",
    "\n",
    "    tt = len(X)\n",
    "    zeros = []\n",
    "    brackets = []\n",
    "\n",
    "    for ii in range(tt - 1):\n",
    "        if f[ii] * f[ii + 1] < 0:  # change of sign\n",
    "            if ii != 0 and ii != tt - 2:  # not at an extreme\n",
    "                if abs(f[ii - 1] - f[ii + 2]) > abs(f[ii] - f[ii + 1]):  # not an asymptote\n",
    "                    brackets += [[X[ii], X[ii + 1]]]\n",
    "            else:\n",
    "                brackets += [[X[ii], X[ii + 1]]]\n",
    "    \n",
    "    sols = []\n",
    "    for br in brackets:\n",
    "        sols.append(optimize.root_scalar(zero_func, args=(V, L),\n",
    "                    bracket=br, method='brentq'))\n",
    "\n",
    "    return [a.root for a in sols]\n",
    "\n",
    "class Mode:\n",
    "    pass\n",
    "\n",
    "\n",
    "def find_modes(a=8.2 / 2, Na=0.12, n_cladding=1.444, w=1.55):\n",
    "    \"\"\"Calculates all the modes in the fiber, and puts them in the list of modes.\n",
    "    It also returns the total number of modes, which is higher than len(modes) because some modes\n",
    "    have degeneracy 2 (L=0) and some have degeneracy 4 (L>0)\"\"\"\n",
    "\n",
    "    k0 = 2 * pi / w\n",
    "    n_core = np.sqrt(Na ** 2 + n_cladding ** 2)\n",
    "\n",
    "    r = np.linspace(0, 3 * a, num=300)\n",
    "    i_rad = round(np.interp(a, r, range(len(r))))  # i-th element at the radius\n",
    "\n",
    "    V = k0 * a * Na\n",
    "    phi = np.linspace(1E-10, pi / 2 - 1E-10, 5000)\n",
    "    X = V * np.sin(phi)\n",
    "    Y = V * np.cos(phi)\n",
    "\n",
    "    solutions = True\n",
    "    L = 0\n",
    "    M = 1\n",
    "    modes = []\n",
    "    tot_modes = 0\n",
    "    while solutions == True:\n",
    "        with np.errstate(invalid='ignore'):\n",
    "            Lhs_1 = X * spe.jv(L + 1, X) / spe.jv(L, X)\n",
    "            Rhs_1 = Y * spe.kv(L + 1, Y) / spe.kv(L, Y)\n",
    "\n",
    "        # sols=CO.FindZerosInterp(Lhs_1-Rhs_1, X)\n",
    "        # my_func = lambda x: (X*spe.jv(L+1,X)/(spe.jv(L,X))-\n",
    "        #                    np.sqrt(V**2-X**2)*spe.kv(L+1,np.sqrt(V**2-X**2))/(spe.kv(L,np.sqrt(V**2-X**2))))\n",
    "        # sols = CO.FindZerosExact(zero_func,x_min = 0, x_max=V, points_grid = 5000)\n",
    "\n",
    "        sols = find_zeros_exact(X, Y, V, L)\n",
    "        for sol in sols:\n",
    "            kt = sol / a\n",
    "            gamma = np.sqrt(V ** 2 - sol ** 2) / a\n",
    "            Er = spe.jv(L, kt * r)\n",
    "            Er[i_rad::] = spe.kv(L, gamma * r[i_rad::]) / spe.kv(L,\n",
    "                    gamma * r[i_rad]) * spe.jv(L, kt * r[i_rad])\n",
    "            Er = Er / max(abs(Er))\n",
    "\n",
    "            mode = Mode()\n",
    "            mode.X = sol\n",
    "            mode.L = L\n",
    "            mode.M = M\n",
    "            mode.Er = Er[:]\n",
    "            mode.V = V\n",
    "            mode.a = a\n",
    "            \n",
    "            mode.neff = np.sqrt(n_core ** 2 - (kt / k0) ** 2)\n",
    "            mode.label = f\"LP({L},{M})\"\n",
    "\n",
    "            mode.degeneracy = 2 if L == 0 else 4\n",
    "\n",
    "            modes.append(mode)\n",
    "\n",
    "            tot_modes += mode.degeneracy\n",
    "\n",
    "            M += 1\n",
    "\n",
    "        M = 1\n",
    "        L += 1\n",
    "        if len(sols) == 0:\n",
    "            solutions = False\n",
    "            \n",
    "    return modes, r, tot_modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-11T23:40:34.455955Z",
     "start_time": "2021-02-11T23:40:33.822952Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3b2762f0f4c94602b38f7f7b7527c1ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "GridspecLayout(children=(Label(value='Core diameter ($\\\\mu$m)', layout=Layout(grid_area='widget001')), Label(v…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1f29ecff216f4e24a9d8d65bf7a8ab81",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Box(layout=Layout(display='flex', flex_flow='column', height='', width='600px'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "slider_diam = widgets.FloatSlider(min=0.1, max=80, value=8.2)\n",
    "slider_Na = widgets.FloatSlider(min=0.0001, max=0.5, step=0.01, value=0.12)\n",
    "slider_lambda = widgets.FloatSlider(min=0.5, max=2.0, step=0.01, value=1.55)\n",
    "\n",
    "text_n_core = widgets.Text()\n",
    "text_results = widgets.Text(value='')\n",
    "btn_calc = widgets.Button(description='Calculate modes')\n",
    "btn_calc.style.button_color = 'lightgray'\n",
    "\n",
    "n_cladding = 1.444\n",
    "\n",
    "def update_text(change):\n",
    "    V = 2 * pi / slider_lambda.value * slider_Na.value * slider_diam.value / 2\n",
    "    \n",
    "    n_core = np.sqrt(slider_Na.value ** 2 + n_cladding ** 2)\n",
    "    \n",
    "    text_n_core.value = f\"n cladding = 1.444, n core = {n_core:.3f}, V = {V:.3f}\"\n",
    "\n",
    "    \n",
    "update_text(0)\n",
    "\n",
    "slider_Na.observe(update_text, names = 'value')    \n",
    "slider_lambda.observe(update_text, names = 'value')\n",
    "slider_diam.observe(update_text, names = 'value')\n",
    "\n",
    "\n",
    "a = slider_diam.value / 2.0\n",
    "Na = slider_Na.value\n",
    "\n",
    "w = slider_lambda.value\n",
    "V = 2 * pi / w * a * Na\n",
    "\n",
    "# diam_core = 50\n",
    "# Na=0.22 #typ. MM fiber\n",
    "    \n",
    "    \n",
    "def fig_mode_profile(mode):\n",
    "\n",
    "    X = mode.X\n",
    "    L = mode.L\n",
    "    V = mode.V\n",
    "    a = mode.a\n",
    "    kt = X / a\n",
    "    gamma = np.sqrt(V ** 2 - X ** 2) / a\n",
    "\n",
    "    E_core = spe.jv(L, kt * a * r_mesh) * np.cos(L * phi_mesh)\n",
    "    E_clad = spe.kv(L, gamma * a * r_mesh) * np.cos(L * phi_mesh) / spe.kv(L, gamma * a) * spe.jv(L, kt * a)\n",
    "    E = E_core * in_core_mesh + E_clad * in_clad_mesh\n",
    "\n",
    "    fig2 = bqplt.figure(fig_margin=dict(top=10, bottom=10, left=10,\n",
    "                        right=10))\n",
    "    fig2.layout.height = '140px'\n",
    "    fig2.layout.width = '140px'\n",
    "    bqplt.heatmap(E, x=xx, y=yy, cmap='RdBu')\n",
    "    line = bqplt.plot(x_core_shape, y_core_shape, 'k', stroke_width=.5)\n",
    "    max_E = np.amax(abs(E))\n",
    "    fig2.axes[0].scale.min = -max_E\n",
    "    fig2.axes[0].scale.max = max_E\n",
    "    for aa in fig2.axes:\n",
    "        aa.visible = False\n",
    "\n",
    "    return fig2\n",
    "\n",
    "\n",
    "def fig_rad_plot(mode, r):\n",
    "\n",
    "    min_y = min(mode.Er)\n",
    "    max_y = max(mode.Er)\n",
    "    xs = bqplot.LinearScale(min=0, max=3 * mode.a)\n",
    "    ys = bqplot.LinearScale(min=min_y, max=max_y)\n",
    "\n",
    "    line1 = bqplt.Lines(x=r, y=mode.Er, scales={'x': xs, 'y': ys})\n",
    "    \n",
    "    grid_line = bqplt.Lines(x=[mode.a, mode.a], y=[min_y, max_y], scales={'x': xs, 'y': ys}, stroke_width=0.5, colors=['black'])\n",
    "    zero_line = bqplt.Lines(x=[0, 3 * mode.a], y=[0, 0], scales={'x': xs, 'y': ys}, stroke_width=0.5, colors=['black'])\n",
    "    \n",
    "    x_ax = bqplt.Axis(orientation='horizontal', scale=xs, tick_values=[], num_ticks=0, visible=False)\n",
    "    y_ax = bqplt.Axis(orientation='vertical', scale=ys, tick_values=[], num_ticks=0)\n",
    "    \n",
    "    new_fig = bqplt.figure(axes=[x_ax, y_ax], marks=[line1, grid_line, zero_line])\n",
    "    new_fig.fig_margin = dict(top=10, bottom=10, left=10, right=10)\n",
    "    new_fig.layout.height = '140px'\n",
    "    new_fig.layout.width = '140px'\n",
    "\n",
    "    # new_fig = bqplt.figure(fig_margin = dict(top=10, bottom=10, left=10, right=10))\n",
    "    # new_fig.layout.height = '100px'\n",
    "    # new_fig.layout.width = '100px'\n",
    "    # bqplt.plot(r,mode.Er)\n",
    "    # bqplt.plot(r,0*r,'k')\n",
    "    # bqplt.plot([mode.a,mode.a],[min(mode.Er),max(mode.Er)],'k:')\n",
    "    # for aa in fig.axes:\n",
    "    #    aa.visible = False\n",
    "\n",
    "    return new_fig\n",
    "\n",
    "\n",
    "def update_table(modes,r):\n",
    "    global num_modes_show\n",
    "    \n",
    "    for ii in range(num_modes_show):\n",
    "        #mode_label.value = modes[0].label\n",
    "        mode = modes[ii]\n",
    "        \n",
    "        show_text = f'Mode {ii+1}:\\n{mode.label}\\nDegeneracy {mode.degeneracy}\\nneff = {mode.neff:0.6f}'\n",
    "        \n",
    "        text_label = widgets.Textarea(value=show_text)\n",
    "        text_label.layout.display = 'flex'\n",
    "        text_label.layout.height = '140px'\n",
    "        \n",
    "        mode_row_list[ii].children = [text_label, fig_rad_plot(modes[ii], r), fig_mode_profile(modes[ii])]\n",
    "        \n",
    "    box_modes.children = mode_row_list[:len(modes)]    \n",
    " \n",
    "    return\n",
    "\n",
    "num_modes_show = 0\n",
    "\n",
    "def btn_eventhandler(obj):\n",
    "    plot()\n",
    "\n",
    "\n",
    "def plot():\n",
    "    global num_modes_show\n",
    "\n",
    "    box_modes.children = []\n",
    "    for ii in range(num_modes_show):\n",
    "        bqplt.close(mode_row_list[ii].children[1])\n",
    "        bqplt.close(mode_row_list[ii].children[2])\n",
    "\n",
    "    w = slider_lambda.value\n",
    "    diam_core = slider_diam.value  # um, typ. SMF28\n",
    "    Na = slider_Na.value  # SMF-28\n",
    "    a = diam_core / 2\n",
    "\n",
    "    (modes, r, tot_modes) = find_modes(a=a, Na=Na, w=w,\n",
    "            n_cladding=n_cladding)\n",
    "    modes.sort(key=lambda x: x.neff, reverse=True)\n",
    "\n",
    "    num_modes = len(modes)\n",
    "    text_results.value = f'Distinct modes: {num_modes}. Total modes: {tot_modes}'\n",
    "    num_modes_show = num_modes\n",
    "    \n",
    "    if num_modes_show > 100:\n",
    "        num_modes_show = 100\n",
    "    \n",
    "    update_table(modes, r)\n",
    "\n",
    "\n",
    "grid = widgets.GridspecLayout(4, 20)\n",
    "\n",
    "grid[0, :4] = widgets.Label('Core diameter ($\\mu$m)')\n",
    "grid[1, 0:4] = widgets.Label('Num. Aperture (NA)')\n",
    "grid[2, :4] = widgets.Label('Wavelength ($\\mu$m)')\n",
    "grid[0, 4:10] = slider_diam\n",
    "grid[1, 4:10] = slider_Na\n",
    "grid[2, 4:10] = slider_lambda\n",
    "grid[3, :4] = btn_calc\n",
    "grid[3, 4:10] = text_results\n",
    "grid[1, 11:] = text_n_core\n",
    "\n",
    "mode_label = widgets.Text()\n",
    "\n",
    "mode_row_layout = widgets.Layout(width='600px', height='150px',\n",
    "                                 border='solid 1px')\n",
    "mode_row_list = []\n",
    "for ii in range(100):\n",
    "    mode_row_list += [widgets.HBox([mode_label, mode_label,\n",
    "                      mode_label], layout=mode_row_layout)]\n",
    "\n",
    "mode_list_layout = widgets.Layout(width='600px', height='',\n",
    "                                  flex_flow='column', display='flex')\n",
    "box_modes = widgets.Box(children=[], layout=mode_list_layout)\n",
    "\n",
    "display(grid)\n",
    "display(box_modes)\n",
    "\n",
    "btn_calc.on_click(btn_eventhandler)\n",
    "plot()"
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}