Skip to content

Instantly share code, notes, and snippets.

@ELC

ELC/environment.yml

Last active February 11, 2021 23:45
Show Gist options
  • Save ELC/4971a4ab3fe64fbdee3875e6c724e58c to your computer and use it in GitHub Desktop.
Save ELC/4971a4ab3fe64fbdee3875e6c724e58c to your computer and use it in GitHub Desktop.
Download ZIP
fiber
Raw
environment.yml
dependencies:
- python=3.7
- ipywidgets==7.5.1
- bqplot==0.11.0
- numpy== 1.19.2
- scipy==1.5.2
Display the source blob
Display the rendered blob
Raw
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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment