fiber_modes.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FIMOC: Fiber-optic Mode Online Calculator\n",
    "\n",
    "\n",
    "Welcome to FIMOC, the Fiber Mode Online Calculator, an online calculator of fiber modes. Here you will be able to calculate a visualize the propagating modes of any step-index fiber of your choice, as in this screenshot:\n",
    "\n",
    "![fig_FIMOC_fiber_modes_screenshot](fiber_modes_images/FIMOC_screenshot.png)\n",
    "\n",
    "*Screenshot of FIMOC showing the modes of a multimode fiber*\n",
    "\n",
    "If you want to go directly to the software, scroll to the bottom, but if you are interested in where these modes come from, read on.\n",
    "\n",
    "## Optical fibers\n",
    "\n",
    "**Did you know that light in optical fibers follow the same patterns as the mechanical vibrations of a drum?** In the next sections, we'll see why. Let's start from the beginning. Optical fibers are long cylinders of transparent glass or plastic, which are characterized by having a refractive index variation along its radial direction. The central region is called the *core* of the fiber, while the external one is called the *cladding*, which is typically made of the same material with a different *doping* concentration. In general, the *core* region has a higher refractive index that the *cladding*, because light hitting an interface from higher to lower index can undergo total internal reflection, which is a necessary condition for guiding, or at least for guiding with a reasonable loss value.\n",
    "\n",
    "A very typical type of fiber is called the *step-index* fiber, which has constant refractive index values at the core and cladding:\n",
    "\n",
    "$$\n",
    "n(r)= \\left\\{\n",
    "\\begin{array}{ll}\n",
    "        n_{core} & r\\le a  \\\\\n",
    "        n_{clad} & r>a\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "where $a$ is the radius of the core, and $n_{core}$ and $n_{clad}$ are the indices of the core and the cladding, respectively. The radius of the cladding doesn't really matter much, because guided light cannot propagate along the cladding, it rapidly decays along the first few microns of cladding, so to simplify the problem, we can consider the cladding radius to be infinite.\n",
    "\n",
    "## Numerical aperture\n",
    "\n",
    "A very important parameter of a fiber is its numerical aperture, which is a number that tells you about its acceptance angle when coupling light from outside. Precisely, it is equal to the sine of the half-angle of the acceptance cone.  If you want the light to be confined in the core, you have to couple the light so that when it propagates in the core within the fiber's critical angle. In a step-index fiber, you can easily calculate (by applying Snell's law) how much this angle should be (within the ray-optics approximation, which is quite valid for multimode fibers), so the numerical aperture of the fiber turns out to be:\n",
    "\n",
    "$$\n",
    "\\text{NA} = \\sqrt{n_{core}^2-n_{clad}^2}\n",
    "$$\n",
    "\n",
    "(you can find a demonstration [here](https://vlab.amrita.edu/?sub=1&brch=189&sim=343&cnt=1) for instance). When you buy a fiber, the datasheet will most likely report the numerical aperture of the fiber, which tells you not only about the collecting (or emitting) angle of the fiber but also about the fiber's index contrast, so remember this definition, as we will need this parameter later on.\n",
    "\n",
    "Another information you will find about the fiber you are considering to buy is whether the fiber is *single-mode* or *multi-mode*. *Modes? What are you talking about??* Bear with me, things will now get more interesting. \n",
    "\n",
    "\n",
    "\n",
    "## Electromagnetic wave equations\n",
    "\n",
    "When light propagates in the fiber, geometrical optics and the *critical angle* concept provides a very useful picture, but when things get small, geometrical optics fails to provide a precise model. If you want the full picture you need to solve the actual electromagnetic wave equations, so that's what we will do now. Beware, there's some math ahead, but don't be scared, you don't really have to understand every single equation, I just want you to have all the details in case you are interested, but feel free to skip this section if you feel overwhelmed.\n",
    "\n",
    "So let's start at the beginning. When light propagates through a medium with a refractive index function $n(\\vec{r})$, any of the components of the electric or magnetic field follows the so-called Helmholtz's equation [1]:\n",
    "\n",
    "$$\n",
    "\\nabla^2 U + n^2(\\vec{r})k_0^2 U=0\n",
    "$$\n",
    "\n",
    "where $n$ is the refractive index and $k_0=2\\pi/\\lambda_0$, being $\\lambda_0$ the wavelength in vacuum.\n",
    " \n",
    "Now let us assume that our medium has cylindrical symmetry. That means that it is convenient to apply Helmholtz's equation in cylindrical coordinates:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial^2U}{\\partial r^2}+\\frac{1}{r}\\frac{\\partial U}{\\partial r}+\\frac{1}{r^2}\\frac{\\partial ^2U}{\\partial\\phi ^2}+\\frac{\\partial ^2U}{\\partial z^2}+n(r)^2k_0^2U=0\n",
    "$$\n",
    "\n",
    "To solve this equation in 3D let's try with a trial solution separating the variables:\n",
    "\n",
    "\n",
    "$$U(r,\\phi,z)=u(r)e^{-jl\\phi}e^{-j\\beta z}$$  for  $$l=0,\\pm 1,\\pm 2...$$\n",
    "\n",
    "where $\\beta$ is the propagating constant along $z$ and the integer number $l$ represents the azimutal dependence of the mode: when $l=0$, it does not depend on the azimut angle $\\phi$, and when $l\\ne0$, you can imagine the 2D mode as a cake subdivided by cuts through its center, so $l$ is simply the number of cuts.\n",
    "\n",
    "When we substitute this differential equation in the previous expression, we get:\n",
    "\n",
    "$$\n",
    "\\frac{d^2u}{dr^2}+\\frac{1}{r}\\frac{du}{dr}+\\left(n^2(r)k_0^2-\\beta^2-\\frac{l^2}{r^2}\\right)u=0 \n",
    "$$\n",
    "\n",
    "which is the differential equation we have to solve for the specific radial index profile of our fiber.\n",
    "\n",
    "The differential equation shown before can have many solutions. Every solution is called a *mode*, which is represented by the 2D profile of the electromagnetic fields that can propagate through the fiber satifying the boundary conditions. Same happens for example with vibration modes of a string under tension, each vibration harmonic is a mode. In our case, for every $l$ parameter, you may find zero, one, or many possible modes, each one characterized by a different propagation constant $\\beta$. Some fibers have a single propagating mode in total, which are called *single-mode* fibers, and some have more than one mode, which are called, what a surprise, *multimode* fibers. \n",
    "\n",
    "\n",
    "\n",
    "## Solution for step-index fibers\n",
    "\n",
    "So the case we are going to solve is the step-index fiber, which we previously defined. The good news is that this equation has analytic solution, which is:\n",
    "\n",
    "$$\n",
    "u(r) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "        J_l(k_T r) & r\\le a  \\\\\n",
    "        K_l(\\gamma r) & r>a\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "where $k_T$ and $\\gamma$ are defined as:\n",
    "\n",
    "$$\n",
    "k_T^2=n_{core}^2k_0^2-\\beta^2 \\\\\n",
    "\\gamma^2=\\beta^2-n_{clad}^2k_0^2\n",
    "$$\n",
    " \n",
    "and J and K are the [Bessel functions](https://en.wikipedia.org/wiki/Bessel_function) of the first and the modified second kind respectively. If it's the first time you hear about Bessel functions don't panic, they are functions that appear very often as solutions to differential equations in cylindrical or spherical coordinates, and many scientific calculators or numerical packages have them in their arsenal. You can simply consider the Bessel functions of the first kind as oscillating functions with a decaying amplitude, and the ones of the modified second kind as some kind of decaying stretched exponential. They look like this:\n",
    "\n",
    "![fig_Bessel_functions](fiber_modes_images/Bessel_functions.png)\n",
    "\n",
    "*Bessel functions of the first kind J and modified second kind K*\n",
    "\n",
    "Actually, the Bessel functions of the first kind are the solutions of the mechanical vibration of a drum with a circular shape, as shown [here](https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane). The only difference with the fiber-optic case is that the decay at the edges is abrupt for the drum case, and is gradual with a modified Bessel function in the fiber-optic case, but in the core are the modes follow the same mathematical function.\n",
    "\n",
    "So you may think that we have found the solution, so we can all go home, but unfortunately we are not done yet. We still need to calculate the propagation constant and the actual electric and magnetic field components that satisfy the boundary conditions, and that is not straightforward in the general case. However in realistic fibers, we can make an aproximation which will simplify the task.\n",
    "\n",
    "## Linearly polarized (LP) modes\n",
    "\n",
    "Most fibers in the real world have a weak refractive index contrast ($n_{core}-n_{clad}<<1$), which is called *weakly guiding* condition. In this approximation, the boundary conditions at the interface can be simplified so that both electric and magnetic fields are continuous and derivable at $r=a$. This greatly simplifies the calculation of the modes, although the solution of the actual propagation constants is not analytical but has to be calculated numerically. Once we get the scalar function $u(r)$, the small index contrast allows us to apply the $u(r)$ profile to any electric field component to get the actual electric field profiles for each polarization, which we choose vertical and horizontal for convencience. Thus in this approximation, each LP mode actually corresponds to two differently polazired modes, one vertical and one horizontal, with identical mode profiles and propagating properties.\n",
    "\n",
    "\n",
    "So this program numerically calculates all possible solutions that satisfy the boundary conditions, and computes the mode profiles for each case. But how many modes does a certain fiber have? To answer that question, it is very helpful to define some new variables (don't worry, we're almost there, I promise). First let us define a normalized adimensional parameter called the $V$-number or the $V$-parameter of the fiber:\n",
    "\n",
    "$$\n",
    "V=2\\pi\\frac{a}{\\lambda_0}\\text{NA}\n",
    "$$\n",
    "\n",
    "This number is very important because it is a single number that can tell you about the number of modes the fiber can hold at a certain wavelength.\n",
    "\n",
    "In the following plot you can see how many modes there are as a function of the $V$ number.\n",
    "\n",
    "![fig_modes_vs_V_number](fiber_modes_images/fiber_modes_vs_V.png)\n",
    "\n",
    "*Calculated modes vs. V-number.*\n",
    "\n",
    "When V < 2.405, the fiber only has one mode (it is a single-mode fiber). Beyond V=2.405, a higher-order mode emerges. That means that V=2.405 is called the *cut-off* of the higher order mode, below which that mode cannot propagate. Every high-order-mode has a *cut-off*, but the first cut-off is the most important because it determines whether the fiber is single-mode or multimode.\n",
    "\n",
    "It is worth noting that the V-number, thus the single- or multi-mode condition, not only depends on the fiber geometrical properties, but also on the wavelength. So a single-mode fiber always becomes multimode at some shorter wavelength.\n",
    "\n",
    "You may notice that the fundamental mode has no cut-off, does it mean that a fiber can guide any wavelength? Mathematically, yes, but in practice no, because some approximations fall apart when guiding becomes extremely weak (V << 1). In particular, when V<<1:\n",
    "\n",
    "- bending loss becomes extremely high (even for very smooth curves)\n",
    "- the radius of the cladding becomes important, as the evanescent field in the cladding can reach very far, making the infinite cladding approximation invalid.\n",
    "- silica glass absorbs wavelengths above $\\gtrsim 2 \\mu$m anyway.\n",
    "\n",
    "\n",
    "Now that we understand the mode cut-offs let's see how LP modes look like. They are characterized by two parameters, the $l$ parameter, as we defined before (the number of cuts on the cake through its center), and a second parameter $m$ which is the number of zeros when going along the radial direction:\n",
    "\n",
    "![fig_LP_modes](fiber_modes_images/fiber_modes_white.png)\n",
    "\n",
    "*Example of some LP modes in a typical multimode fiber. Red and blue colors represent electric fields with opposite signs.*\n",
    "\n",
    "\n",
    "## Number of modes\n",
    "\n",
    "It is worth mentioning that every LP mode can actually represent more than one propagating mode. Let's remember that each mode can have two polarizations (like the vibrations of a string, by the way), which in the low-index-contrast approximation they can be considered as degenerate. But when $|l|$>0, the mode also has two helical polarities, corresponding to $+l$ and $-l$. This means that modes with $l$=0 have a degeneracy of 2 while modes with $|l|$>0 have degeneracy 4. However they are not typically represented as they have the same shape. So when speaking about the number of modes in a fiber, it is always a good idea to clarify if you are talking about distinct LP modes or total number of modes including degenerates.\n",
    "\n",
    "A useful equation to estimate the total number of modes (including degenerates) as a function of $V$ is given by this equation (demonstrated in [Saleh's book](#Bibliography)):\n",
    "\n",
    "$$\n",
    "N_{modes} \\approx \\frac{V^2}{2} \\text{    ( for $V$ >> 1)} \n",
    "$$\n",
    "\n",
    "which becomes precise for high values of $V$ (typical multimode fibers).\n",
    "\n",
    "\n",
    "## Modal dispersion\n",
    "\n",
    "If you want to transmit data through a multimode fiber, there is a problem: each mode not only has a different 2D profile, but it also propagates at a different velocity. This is because the velocity of each mode is related to the propagation constant $\\beta$ of that specific mode, so it is different for every mode. Actually, the propagation constant by itself just tells you the phase velocity of the mode, which does not coincide with the actual mode propagation velocity which is called group velocity. This number also depends on how the phase velocity depends on wavelength. It can be calculated with these equations:\n",
    "\n",
    "\n",
    "$$\n",
    "v_{g} = \\frac{c}{n_g}\n",
    "$$\n",
    "\n",
    "where $v_{g}$ is the group velocity, and $n_g$ is the group index, defined as:\n",
    "\n",
    "$$\n",
    "n_g = n_{eff} - \\lambda \\frac{dn_{eff}}{d\\lambda}\n",
    "$$\n",
    "\n",
    "FIMOC also calculates the group velocity of each mode. It does so by calculating the mode spectrum not only in the required wavelength, but also at slightly shifted wavelengths. For doing so it also requires the group index of the material itself, which is given as a parameter. The program does the approximation of considering the same index vs. wavelength slope on the cladding and on the core, to avoid requiring too many parameters.\n",
    "\n",
    "It is worth noting that although the fundamental mode has the highest effective index of all modes (the one closest to the core index), it has the lowest group index of all modes (with the exception of very weakly-guided modes that may appear in some cases). This is reasonable as in geometrical optics, higher order modes must propagate slower that the fundamental mode as they travel a longer distance because they rebound many times at the interfaces.\n",
    "\n",
    "Considering the range of group velocities within the modal spectrum, it is expected that a short pulse of light will spread out after a certain length due to this speed variability, which is called modal dispersion. When transmitting information, each bit can spread and intermix with the adjacent ones up to the point of making the signal unreadable. For this reason, modal dispersion imposes a bitrate limitation in multimode fibers. With this program you can easily calculate it by taking the group velocity range of all the modes.\n",
    "\n",
    "This means that if you want to transmit a high-bitrate signal though a long distance, you should use a single-mode fiber for this reason. I'm not saying that single mode fibers have no dispersion at all, there is also some dispersion in those, which is called chromatic dispersion, but that's a different story.\n",
    "\n",
    "Another interesting fact is that you can greatly reduce modal dispersion in a multimode fiber by gradually doping the refractive index profile with a parabolic shape, rather than a step-index profile, but that's also a different story. Let's get back to step-index fibers, which are the ones that this software deals with.\n",
    "\n",
    "\n",
    "\n",
    "## FIMOC interface\n",
    "\n",
    "FIMOC calculates the modes of a step-index fiber, given its core diameter, numerical aperture and wavelength. These parameters can be set using the sliders. When the sliders are moved, the V-parameter of the fiber is also shown in real time, to have an initial idea of how many modes you can expect.\n",
    "\n",
    "Even though it is not used for the mode calculation, the refractive index of the cladding can also be provided by the user. This parameter is only used to calculate the core index (from the fiber NA parameter) and to show the effective index of each mode, defined as $n_{eff} = \\beta/k_0$. The effective index is more meaningful than $\\beta$ because it is a value that is between $n_{core}$ and $n_{clad}$, and tells you that the mode propagates with a phase velocity equivalent to a medium with a refractive index equal to $n_{eff}$. So when $n_{eff}$ is very close to $n_{clad}$, the mode is weakly confined, while a mode with $n_{eff}$ close to $n_{core}$ is a mode that comfortably fits in the core with almost no evanescent field in the cladding. The program also takes the group velocity of the cladding to calculate the group velocity of each mode, assuming that the core material has the same shift between refractive index and group index.\n",
    "\n",
    "You can play with the parameters and see the modes for any case you want. At the moment there is a restriction on the highest V number to 30, as beyond that value the number of modes becomes more that 100, which may take long time to plot. Besides, when your fiber can hold more than 100 modes, you don't really need to model light propagation as a combination of modes, but the ray-optics model is probably good enough in most cases.\n",
    "\n",
    "You may want to try some typical fiber parameters:\n",
    "- **Standard single-mode fiber SMF-28** (single-mode in range $\\lambda$ > 1.25$\\mu$m):\n",
    "    - Core diameter: 8.2 $\\mu$m\n",
    "    - NA = 0.12\n",
    "    - $\\lambda$ = 1.55 $\\mu$m\n",
    "- **A typical multimode fiber**:\n",
    "    - Core diameter: 50 $\\mu$m\n",
    "    - NA = 0.22\n",
    "    - $\\lambda$ = 1.55 $\\mu$m\n",
    "    \n",
    "Click any of the widget buttons below to launch FIMOC. (Sometimes it may take long for the server to launch, sometimes even more than 20-30 s. so please be patient...)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-11T23:40:34.455955Z",
     "start_time": "2021-02-11T23:40:33.822952Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "da19682e57754dd89c4aed8f07e20062",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(VBox(children=(Label(value='Core diameter (um)', layout=Layout(display='flex', justify_content=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0b327fa674dd48daa2ee1cc011a23d42",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3f024cc198d940228cd28bc5bd7b94d5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Button(description='PLOT EVERY MODE', layout=Layout(visibility='hidden'), style=ButtonStyle(but…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "792876e9f07b41c18f4b1c709cae326f",
       "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": [
    "## Imports\n",
    "\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import scipy.special as spe\n",
    "from scipy import optimize\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import HBox, VBox, Layout, Label, IntSlider, FloatSlider, Button, HTML, Output, FloatText, Box\n",
    "\n",
    "%matplotlib inline\n",
    "plt.style.use(\"bmh\")\n",
    "\n",
    "pi = np.pi\n",
    "\n",
    "## Initialization\n",
    "dpi = 100\n",
    "plt.rcParams[\"figure.dpi\"] = dpi\n",
    "\n",
    "mesh_limit = 2.\n",
    "xx = np.linspace(-mesh_limit, mesh_limit, 60)\n",
    "yy = np.linspace(-mesh_limit, mesh_limit, 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",
    "\n",
    "in_core_mesh = ones_mesh.copy()\n",
    "in_core_mesh[r_mesh > 1] = zeros_mesh[r_mesh > 1]  # mask core with ones\n",
    "\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 perimeter 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\n",
    "\n",
    "## Functions\n",
    "\n",
    "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",
    "        \n",
    "        if f[ii] * f[ii + 1] < 0:  # change of sign\n",
    "            \n",
    "            if ii != 0 and ii != tt - 2:  # not at an extreme\n",
    "                \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",
    "                    \n",
    "            else:\n",
    "                brackets += [[X[ii], X[ii + 1]]]\n",
    "    \n",
    "    sols = []\n",
    "    \n",
    "    for br in brackets:\n",
    "        optimum = optimize.root_scalar(zero_func, args=(V, L), bracket=br, method='brentq')\n",
    "        sols.append(optimum)\n",
    "\n",
    "    return [a.root for a in sols]\n",
    "\n",
    "\n",
    "def calc_Er(mode,r):\n",
    "    kt = mode.X / mode.a\n",
    "    gamma = np.sqrt(mode.V ** 2 - mode.X ** 2) / mode.a\n",
    "    Er = np.empty([len(r)])\n",
    "    \n",
    "    correction_clad = spe.jv(mode.L, kt * mode.a) / spe.kv(mode.L, gamma * mode.a)\n",
    "    \n",
    "    for ii in range(len(r)):\n",
    "        \n",
    "        if r[ii] < - mode.a:\n",
    "            Er[ii] = spe.kv(mode.L, gamma * abs(r[ii])) * correction_clad\n",
    "            if mode.L % 2:\n",
    "                Er[ii] = -Er[ii]\n",
    "                \n",
    "        elif r[ii] > mode.a:\n",
    "            Er[ii] = spe.kv(mode.L, gamma * r[ii]) * correction_clad\n",
    "            \n",
    "        else:\n",
    "            Er[ii] = spe.jv(mode.L, kt * r[ii])\n",
    "    \n",
    "    Er = Er / max(abs(Er))\n",
    "    return Er\n",
    "\n",
    "def calc_E_mesh(mode,r_mesh,phi_mesh,in_core_mesh,in_clad_mesh):\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",
    "    E = E / np.amax(E)\n",
    "    return E\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",
    "    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",
    "    L = 0\n",
    "    M = 1\n",
    "    modes = []\n",
    "    tot_modes = 0\n",
    "    \n",
    "    with np.errstate(invalid='ignore'):\n",
    "        while True:\n",
    "            sols = find_zeros_exact(X, Y, V, L)\n",
    "        \n",
    "            for sol in sols:\n",
    "                mode = Mode()\n",
    "                mode.X, mode.L, mode.M, mode.V, mode.a = sol, L, M, V, a\n",
    "\n",
    "                mode.neff = np.sqrt(n_core**2 - (sol/(a*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",
    "\n",
    "            if len(sols) == 0:\n",
    "                break\n",
    "            \n",
    "    return modes, tot_modes\n",
    "\n",
    "\n",
    "## Interactive elements\n",
    "\n",
    "## Sliders\n",
    "slider_diam = FloatSlider(min=2.0, max=80, value=50)\n",
    "slider_Na = FloatSlider(min=0.02, max=0.5, step=0.001, value=0.12)\n",
    "slider_lambda = FloatSlider(min=0.5, max=2.0, step=0.01, value=1.55)\n",
    "text_n_clad = FloatText(min = 1.0, max = 2.0, step = 0.01, value = 1.444, layout=Layout(width=\"80px\"))\n",
    "text_ng_clad = FloatText(min = 1.0, max = 2.0, step = 0.01, value = 1.4626, layout=Layout(width=\"80px\"))\n",
    "\n",
    "\n",
    "## Labels\n",
    "label_n_core = Label()\n",
    "\n",
    "label_V = Label()\n",
    "\n",
    "label_V_max = Label(value = r'\\(\\color{red} {V max = 30!}\\)')\n",
    "label_V_max.layout.visibility = 'Hidden'\n",
    "\n",
    "label_text_modes_found = Label()\n",
    "\n",
    "## Buttons\n",
    "btn_calc = Button(description='CALCULATE MODES')\n",
    "btn_calc.style.button_color = 'lightgray'\n",
    "\n",
    "btn_plot = Button(description='PLOT EVERY MODE')\n",
    "btn_plot.style.button_color = 'lightgray'\n",
    "btn_plot.layout.visibility = 'Hidden'\n",
    "\n",
    "btn_plot.label_gen_plots = Label(value = 'Generating plots...')\n",
    "btn_plot.label_gen_plots.layout.visibility = 'Hidden'\n",
    "\n",
    "\n",
    "def update_text(obj):\n",
    "    V = 2 * pi / slider_lambda.value * slider_Na.value * slider_diam.value / 2\n",
    "    if V > 30:\n",
    "        label_V_max.layout.visibility = 'Visible'\n",
    "        obj['owner'].value = obj['old']\n",
    "        V = 2 * pi / slider_lambda.value * slider_Na.value * slider_diam.value / 2\n",
    "\n",
    "\n",
    "    n_core = np.sqrt(slider_Na.value ** 2 + text_n_clad.value ** 2)\n",
    "    label_n_core.value = f\"{n_core:.4f}\"\n",
    "    label_V.value = f\"{V:.4f}\"\n",
    "    btn_plot.layout.visibility = 'Hidden'\n",
    "\n",
    "\n",
    "# Initialize\n",
    "update_text(0)\n",
    "\n",
    "## Observe Changes\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",
    "text_n_clad.observe(update_text, names = 'value')\n",
    "text_ng_clad.observe(update_text, names = 'value')\n",
    "\n",
    "## Update Values\n",
    "a = slider_diam.value / 2.0\n",
    "Na = slider_Na.value\n",
    "w = slider_lambda.value\n",
    "V = 2 * pi / w * a * Na\n",
    "n_clad = text_n_clad.value\n",
    "\n",
    "\n",
    "# diam_core = 50\n",
    "# Na=0.22 #typ. MM fiber\n",
    "\n",
    "    \n",
    "\n",
    "\n",
    "num_modes_show = 0\n",
    "\n",
    "def btn_calc_eventhandler(obj):\n",
    "\n",
    "    label_V_max.layout.visibility = 'Hidden'\n",
    "    label_text_modes_found.value = 'Calculating...'\n",
    "\n",
    "    w = slider_lambda.value\n",
    "    diam_core = slider_diam.value\n",
    "    Na = slider_Na.value\n",
    "    n_cladding = text_n_clad.value\n",
    "    ng_cladding = text_ng_clad.value\n",
    "    n_core = np.sqrt(Na ** 2 + n_cladding ** 2)\n",
    "    ng_core = n_core + (ng_cladding-n_cladding)\n",
    "    delta_n = n_core-n_cladding\n",
    "    a = diam_core / 2\n",
    "    r = np.linspace(-3*a, 3*a, 100)\n",
    "    \n",
    "\n",
    "    (modes, tot_modes) = find_modes(a=a, Na=Na, w=w, n_cladding=n_cladding)\n",
    "    \n",
    "    delta_w = 1E-4\n",
    "    delta_n_clad = (n_cladding-ng_cladding)/w*delta_w\n",
    "    (modes_red, tot_modes_red) = find_modes(a=a, Na=Na, w=w+delta_w, n_cladding=n_cladding+delta_n_clad)\n",
    "    (modes_blue, tot_modes_blue) = find_modes(a=a, Na=Na, w=w-delta_w, n_cladding=n_cladding-delta_n_clad)\n",
    "    \n",
    "    for ii,m in enumerate(modes):\n",
    "        new_ng = m.neff-w*(modes_red[ii].neff-modes_blue[ii].neff)/(2*delta_w)\n",
    "        m.ng = new_ng\n",
    "    \n",
    "\n",
    "    M_max = max(m.M for m in modes)\n",
    "    M_list = np.arange(M_max)+1\n",
    "    L_max = max(m.L for m in modes)\n",
    "    modes.sort(key=lambda x: x.neff, reverse=True)\n",
    "    \n",
    "    label_text_modes_found.value = f'Distinct modes found: {len(modes)}. Total modes: {tot_modes}'\n",
    "    obj.modes, obj.r = modes, r\n",
    "    btn_plot.layout.visibility = 'visible'\n",
    "    obj.mode_range_plot.clear_output(wait=True)\n",
    "    with obj.mode_range_plot:\n",
    "        fig, (ax1,ax2) = plt.subplots(1,2,figsize=(600 / dpi, 250 / dpi))\n",
    "        for M in M_list:\n",
    "            modes_subset= [mode_sub for mode_sub in modes if mode_sub.M==M]\n",
    "            y=[m.neff for m in modes_subset]\n",
    "            #x= L*np.ones(len(y))\n",
    "            x=[m.L for m in modes_subset]\n",
    "            ax1.plot(x,y,'.')\n",
    "            y2=[m.ng for m in modes_subset]\n",
    "            ax2.plot(x,y2,'.')\n",
    "        \n",
    "        \n",
    "     \n",
    "        \n",
    "       \n",
    "        ax1.set_ylim([n_cladding-delta_n*.15, n_core+delta_n*.15])\n",
    "        ax1.axhline(n_cladding,linewidth=1, color='black', linestyle = 'dashed')\n",
    "        ax1.axhline(n_core,linewidth=1, color='black',linestyle = 'dashed')\n",
    "        ax1.annotate('n clad',xy = (0.03,0.04),xycoords='axes fraction', fontsize = 6)\n",
    "        ax1.annotate('n core',xy = (0.03,0.91),xycoords='axes fraction', fontsize = 6)\n",
    "        ax1.set_ylabel('$n_{eff}$',fontsize=15)\n",
    "        ax1.set_xlabel('$L$')\n",
    "        \n",
    "        from matplotlib.ticker import MaxNLocator\n",
    "        ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "        if L_max == 0:\n",
    "            ax1.set_xticks([0.0])\n",
    "        \n",
    "        ax2.axhline(ng_core,linewidth=1, color='black',linestyle = 'dashed')\n",
    "      \n",
    "        trans = ax2.get_yaxis_transform()\n",
    "   \n",
    "        ylims = ax2.get_ylim()\n",
    "        ax2.text(0.8,ng_core+0.05*(ylims[1]-ylims[0]),'ng core', transform = trans, fontsize = 6,verticalalignment = 'bottom')\n",
    "        ax2.set_ylabel('$n_{g}$',fontsize=15)\n",
    "        ax2.set_xlabel('$L$')\n",
    "        ax2.xaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "        if L_max == 0:\n",
    "            ax2.set_xticks([0.0])\n",
    "        \n",
    "        fig.tight_layout()\n",
    "            \n",
    "        plt.show(fig)\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "def btn_plot_eventhandler(obj):\n",
    "    obj.label_gen_plots.layout.visibility = 'Visible'\n",
    "    \n",
    "    for ii in range(len(obj.mode_row_list)-1):\n",
    "        obj.mode_row_list[ii + 1].children[1].clear_output(wait=True)\n",
    "        obj.mode_fig_list[ii].clf()\n",
    "        \n",
    "    fig_titles_layout = Layout(width = '590px', border='None', margin = '10px 10px 10px 10px', justify_content = 'space-around')\n",
    "    labels = [\n",
    "        Label('  LP(L,M) mode description', layout=Layout(width='250px')),\n",
    "        Label('Radial plot', layout=Layout(width='190px')),\n",
    "        Label('2D image plot', layout=Layout(width='150px'))\n",
    "    ]\n",
    "    \n",
    "    obj.mode_row_list = [HBox(labels, layout=fig_titles_layout)]\n",
    "    \n",
    "    obj.mode_fig_list = []\n",
    "    \n",
    "    modes, r = btn_calc.modes, btn_calc.r\n",
    "\n",
    "    for ii in range(len(modes)):\n",
    "        mode = modes[ii]\n",
    "        \n",
    "        show_text = f'<b>Mode {ii+1}</b>:<br>{mode.label}<br>Degeneracy {mode.degeneracy}' \\\n",
    "                    f'<br>neff = {mode.neff:0.6f}<br>ng = {mode.ng:0.6f}'\n",
    "        \n",
    "        text_label = HTML(value=show_text, layout=Layout(width='100px'))\n",
    "        \n",
    "        mode_row = HBox([text_label,Output()],layout=mode_row_layout)\n",
    "        obj.mode_row_list.append(mode_row)\n",
    "        \n",
    "        with mode_row.children[1]:\n",
    "            fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(400 / dpi, 190 / dpi))\n",
    "            \n",
    "            ax1.plot(r, calc_Er(modes[ii], r))\n",
    "            \n",
    "            ax1.set_xticks([-mode.a, mode.a])\n",
    "            ax1.set_xticklabels(['-a','a'])\n",
    "            ax1.set_xlim([r[0],r[-1]])\n",
    "            ax1.set_ylabel('E(r)')\n",
    "            ax1.set_yticks([0])\n",
    "            ax1.set_yticklabels([])\n",
    "            ax1.set_position([0.1,0.1,.4,.8])\n",
    "            \n",
    "            E = calc_E_mesh(mode, r_mesh, phi_mesh, in_core_mesh, in_clad_mesh)\n",
    "            \n",
    "            ax2.imshow(E, cmap='RdBu', vmin = -1., vmax = 1.,\n",
    "                       extent=(xx[0], xx[-1], yy[0], yy[-1]),\n",
    "                       aspect='equal')\n",
    "            \n",
    "            ax2.plot(x_core_shape,y_core_shape,'black',linewidth=0.5)\n",
    "            ax2.set_xticks([])\n",
    "            ax2.set_yticks([])\n",
    "            ax2.set_position([0.6, 0.1, .4, .8])\n",
    "            \n",
    "            plt.show(fig)\n",
    "            obj.mode_fig_list.append(fig)\n",
    "            \n",
    "    obj.box_modes.children = obj.mode_row_list\n",
    "    obj.label_gen_plots.layout.visibility = 'hidden'\n",
    "\n",
    "\n",
    "gr_layout= Layout(width='800px', border='solid 1px',\n",
    "                  grid_template_columns='20% 40% 20% 20%',\n",
    "                  grid_tempalte_rows='20% 20% 20% 20%',\n",
    "                )\n",
    "\n",
    "empty = Label()\n",
    "\n",
    "label_layout = Layout(display='flex', justify_content='flex-end')\n",
    "\n",
    "labels_units = [Label('Core diameter (um)', layout=label_layout), Label('Num. Aperture (NA)', layout=label_layout),\n",
    "                Label('Wavelength (um)', layout=label_layout), btn_calc]\n",
    "\n",
    "sliders = [slider_diam, slider_Na, slider_lambda, \n",
    "           label_text_modes_found]\n",
    "\n",
    "labels_variables = [Label('neff cladding', layout=label_layout),\n",
    "                    Label('ng cladding', layout = label_layout),\n",
    "                    Label('n core', layout=label_layout),\n",
    "                    Label('V number', layout=label_layout)]\n",
    "\n",
    "elements = [text_n_clad, text_ng_clad, label_n_core, HBox([label_V, label_V_max]), empty]\n",
    "\n",
    "columns = [VBox(labels_units), VBox(sliders), VBox(labels_variables), VBox(elements)]\n",
    "\n",
    "grid_ctrl = HBox(columns)\n",
    "\n",
    "btn_calc.mode_range_plot = Output()\n",
    "\n",
    "btn_plot_row = HBox([btn_plot, btn_plot.label_gen_plots])\n",
    "\n",
    "mode_row_layout = Layout(width='590px', border='solid 1px', margin = '3px 3px 3px 3px', justify_content = 'space-around')\n",
    "\n",
    "btn_plot.mode_row_list = []\n",
    "\n",
    "mode_list_layout = Layout(width='600px', height='', flex_flow='column', display='flex')\n",
    "btn_plot.box_modes = Box(children=btn_plot.mode_row_list, layout=mode_list_layout)\n",
    "\n",
    "display(grid_ctrl)\n",
    "display(btn_calc.mode_range_plot)\n",
    "display(btn_plot_row)\n",
    "display(btn_plot.box_modes)\n",
    "\n",
    "btn_calc.modes = []\n",
    "btn_calc.r=0.\n",
    "\n",
    "btn_calc.on_click(btn_calc_eventhandler)\n",
    "btn_plot.on_click(btn_plot_eventhandler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Version\n",
    "\n",
    "This is FIMOC version 2.0, published in March 2022.\n",
    "\n",
    "## About the author\n",
    "\n",
    "I am Claudio Oton, I am a researcher at [Scuola Superiore Sant'Anna](https://www.santannapisa.it) in Pisa, Italy, where I work doing research on fiber-optics. Please let me know if there is any issue in the software, a mistake in the text, or any other suggestion you may have. You may contact me at c.oton{at}santannapisa.it\n",
    "\n",
    "\n",
    "## Acknowledgments\n",
    "\n",
    "Special thanks to Marco Vannucci for introducing me to Jupyter Notebooks, which were used for making the software, and to [Ezequiel Castaño](https://elc.github.io/) for his [tutorial](https://elc.github.io/posts/embed-interactive-notebooks/) to make interactive online programs from Jupyter notebooks, and his personal help with FIMOC.\n",
    "\n",
    "\n",
    "## Bibliography\n",
    "\n",
    "1. \"Fundamentals of Photonics\" by Saleh and Teich (Ed. Wiley). This was my main source for making this software. In this book you will find more details and further references.\n"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}