PLaSK library

Loading...
Searching...
No Matches
fiber.py
Go to the documentation of this file.
1# This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
2# Copyright (c) 2022 Lodz University of Technology
3#
4# This program is free software: you can redistribute it and/or modify
5# it under the terms of the GNU General Public License as published by
6# the Free Software Foundation, version 3.
7#
8# This program is distributed in the hope that it will be useful,
9# but WITHOUT ANY WARRANTY; without even the implied warranty of
10# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11# GNU General Public License for more details.
12
13 import numpy as np
14 import matplotlib.pyplot as plt
15
16 from scipy.special import jv, kv
17 from scipy.optimize import fsolve
18
19
20 class Analytic:
21
22 def __init__(self, R=1., nr=3.5, m=1, lam=1000.):
23 self.R = R
24 self.nr = nr
25 self.m = m
26 self.J = lambda z: jv(m, z)
27 self.K = lambda z: kv(m, z)
28 self.D_J = lambda z: 0.5 * (jv(m-1,z) - jv(m+1, z))
29 self.D_K = lambda z: - kv(m-1, z) - m/z * kv(m, z)
30 self.k0 = 2e3*np.pi / lam
31 self.V = R * self.k0 * (nr**2 - 1)**0.5
32
33 nn = np.linspace(1.00, nr, 4001)[1:-1]
34 dets = np.abs(self.fun fun(nn))
35 self._minima = nn[1:-1][(dets[1:-1] <= dets[:-2]) & (dets[1:-1] < dets[2:])][::-1]
36
37 def fun(self, neff):
38 """
39 Characteristic function for finding fiber modes
40
41 Taken from Snyder & Love, p.253
42 """
43 kz = neff * self.k0
44 U = self.R * ((self.k0 * self.nr)**2 - kz**2)**0.5
45 W = self.R * (kz**2 - self.k0**2)**0.5
46 fcore = self.D_J(U) / (U * self.J(U))
47 fclad = self.D_K(W) / (W * self.K(W))
48 lhs = (fcore + fclad) * (fcore + fclad / self.nr**2)
49 rhs = (self.m * kz / (self.k0 * self.nr))**2 * (self.V / (U * W))**4
50 return lhs - rhs
51
52 def field(self, neff, r):
53 kz = neff * self.k0
54 sr = np.sign(r)
55 r = np.abs(r)
56 m = self.m
57 U = self.R * ((self.k0 * self.nr)**2 - kz**2)**0.5
58 W = self.R * (kz**2 - self.k0**2)**0.5
59 Δ = 0.5 * (self.nr**2 - 1) / self.nr**2
60 JU = jv(m, U)
61 JUp = jv(m+1, U)
62 JUm = jv(m-1, U)
63 KW = kv(m, W)
64 KWp = kv(m+1, W)
65 KWm = kv(m-1, W)
66 if m != 0:
67 b1 = 0.5 / U * (JUm - JUp) / JU
68 b2 = - 0.5 / W * (KWm + KWp) / KW
69 F2 = (self.V / (U*W))**2 * m / (b1 + b2)
70 a1 = 0.5 * (F2 - 1)
71 a2 = 0.5 * (F2 + 1)
72 Er = - np.where(r <= self.R,
73 (a1 * jv(m-1,U*r) + a2 * jv(m+1,U*r)) / JU,
74 U/W * (a1 * kv(m-1,W*r) - a2 * kv(m+1,W*r)) / KW)
75 Ep = + np.where(r <= self.R,
76 (a1 * jv(m-1,U*r) - a2 * jv(m+1,U*r)) / JU,
77 U/W * (a1 * kv(m-1,W*r) + a2 * kv(m+1,W*r)) / KW)
78 Ez = - U / (kz*self.R) * np.where(r <= self.R,
79 jv(m, U*r) / JU,
80 kv(m, W*r) / KW)
81 if m % 2 == 0:
82 Er[sr < 0] = -Er[sr < 0]
83 Ep[sr < 0] = -Ep[sr < 0]
84 else:
85 Ez[sr < 0] = -Ez[sr < 0]
86 else:
87 fcore = self.D_J(U) / (U * JU)
88 fclad = self.D_K(W) / (W * KW)
89 if abs(fcore + fclad) < 1e-6:
90 Ep = - np.where(r <= R, jv(1, U*r) / JUp, kv(1, W*r) / KWp)
91 Er = Ez = np.zeros(len(Ep))
92 else:
93 Er = np.where(r <= R, jv(1, U*r) / JUp, nr**2 * kv(1, W*r) / KWp)
94 Ep = np.zeros(len(Er))
95 Ez = np.where(r <= R, U / (kz*R) * jv(0, U*r) / JUp, - W / (kz*R) * nr**2 * kv(0, W*r) / KWp)
96 M = np.real(Er * np.conj(Er) + Ep * np.conj(Ep) + Ez * np.conj(Ez))
97 fm = 1 / np.max(M)**0.5
98 Er *= fm
99 Ep *= fm
100 Ez *= fm
101 return np.array([Ep, Er, Ez]).T
102
103 def __len__(self):
104 return len(self._minima)
105
106 def __getitem__(self, num):
107 if isinstance(num, slice):
108 return [fsolve(self.fun fun, x)[0] for x in self._minima[num]]
109 else:
110 return fsolve(self.fun fun, self._minima[num])[0]
111
112
113 def make_radial_plot(Er, r=None, m=1, R=1., neff=None, c=1, axs=None):
114 """Plot radial profile
115
116 Args:
117 Er (array of Er): Field to plot.
118 r (array, optional): Radial positions. If None, it is take from data mesh.
119 m (int): Angular mode param.
120 R (float): Fiber radius.
121 """
122 if r is None:
123 r = Er.mesh.axis0
124 Er = Er.array[:,0,:]
125 if axs is None:
126 fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
127 for ax in ax1, ax2:
128 ax.axvline(-R, color='k', ls='--', lw=0.5)
129 ax.axvline(+R, color='k', ls='--', lw=0.5)
130 ax1.plot((np.nan,), (np.nan,), color='k', ls='-', label="$E_r$")
131 ax1.plot((np.nan,), (np.nan,), color='k', ls='-', alpha=0.35, label=r"$E_\varphi$")
132 ax1.plot((np.nan,), (np.nan,), color='k', ls='--', label="$E_z$")
133 ax2.set_xlabel("$r$ (µm)")
134 leg1 = ax1.legend(loc='best', prop={'size': 6})
135 try: leg1.set_draggable(True)
136 except AttributeError: pass
137 fig.set_tight_layout({'pad': 0.1})
138 else:
139 ax1, ax2 = axs
140 ax1.plot(r, Er[:,1], color=f'C{c}', ls='-')
141 ax1.plot(r, Er[:,0], color=f'C{c}', ls='-', alpha=0.5)
142 ax1.plot(r, Er[:,2], color=f'C{c}', ls='--')
143 ax2.plot(r, np.sum(np.real(Er * np.conj(Er)), 1), color=f'C{c}', label=f"{c}: {neff.real:.4f}" if neff is not None else None)
144 leg2 = ax2.legend(loc='best', prop={'size': 6})
145 try: leg2.set_draggable(True)
146 except AttributeError: pass
147 return ax1, ax2
148
149
150 def make_polar_plot(Er, r=None, m=1, R=1.):
151 """Plot polar maps
152
153 Args:
154 Er (array of Er): Field to plot.
155 r (array, optional): Radial positions. If None, it is take from data mesh.
156 m (int): Angular mode param.
157 R (float): Fiber radius.
158 """
159 if r is None:
160 r = Er.mesh.axis0
161 Er = Er.array[:,0,:2]
162
163 p = np.linspace(0., 2*np.pi, 37)
164 rr, pp = np.meshgrid(r, p)
165 phi = np.linspace(0., 2*np.pi, 1001)
166
167 E = Er[None,:,:]
168 E = np.repeat(E, len(p), 0)
169
170 mer = np.max(np.abs(E.ravel().real))
171 mei = np.max(np.abs(E.ravel().imag))
172 if mer >= mei: E = E.real / mer
173 else: E = E.imag / mei
174
175 if m != 0:
176 E[:,:,0] *= np.sin(m * pp)
177 E[:,:,1] *= np.cos(m * pp)
178
179 F = np.sum((E * E.conj()).real, 2)
180
181 fs = plt.rcParams['figure.figsize'][0]
182 fig, ax = plt.subplots(subplot_kw=dict(projection='polar'), figsize=(fs,fs))
183 ax.tick_params('y', colors='w')
184 ax.tick_params(axis='x', left='off', top='off', right='off', bottom='off', labelleft='off', labeltop='off', labelright='off', labelbottom='off')
185 ax.tick_params(axis='y', left='off', top='off', right='off', bottom='off', labelleft='off', labeltop='off', labelright='off', labelbottom='off')
186 ax.set_xticklabels([])
187 ax.set_yticklabels([])
188 plt.grid(False)
189
190 plt.contourf(pp, rr, F, 64)
191
192 Ex = E[:,:,1] * np.cos(pp) - E[:,:,0] * np.sin(pp)
193 Ey = E[:,:,1] * np.sin(pp) + E[:,:,0] * np.cos(pp)
194 plt.quiver(pp, rr, Ex, Ey, pivot='mid', scale=30, color='w')
195
196 plt.plot(phi, (R,)*len(phi), 'w', lw=0.5)