PLaSK library
Loading...
Searching...
No Matches
vcsel.py
Go to the documentation of this file.
1#!/usr/bin/env python3
2# This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3# Copyright (c) 2022 Lodz University of Technology
4#
5# This program is free software: you can redistribute it and/or modify
6# it under the terms of the GNU General Public License as published by
7# the Free Software Foundation, version 3.
8#
9# This program is distributed in the hope that it will be useful,
10# but WITHOUT ANY WARRANTY; without even the implied warranty of
11# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12# GNU General Public License for more details.
13
14import unittest
15
16import sys
17
18from numpy import *
19
20from plask import *
21from plask import material, geometry, mesh
22from optical import modal
23
24
25plask.config.axes = "rz"
26
27
28class VCSEL(unittest.TestCase):
29
30 def setUp(self):
31 plask.config.axes = 'rz'
32 self.manager = plask.Manager()
33 self.manager.load('''
34 <plask loglevel="debug">
35 <materials>
36 <material name="GaAs" base="semiconductor">
37 <nr>3.53</nr>
38 <absp>0.</absp>
39 </material>
40 <material name="AlGaAs" base="semiconductor">
41 <nr>3.08</nr>
42 <absp>0.</absp>
43 </material>
44 <material name="AlAs" base="semiconductor">
45 <nr>2.95</nr>
46 <absp>0.</absp>
47 </material>
48 <material name="AlOx" base="semiconductor">
49 <nr>1.53</nr>
50 <absp>0.</absp>
51 </material>
52 <material name="InGaAs" base="semiconductor">
53 <nr>3.53</nr>
54 <absp>0.</absp>
55 </material>
56 </materials>
57 <geometry>
58 <cylindrical axes="rz" name="vcsel" outer="extend" bottom="GaAs">
59 <stack name="layers">
60 <block dr="5." dz="0.06949" material="GaAs"/>
61 <stack name="top-dbr" repeat="24">
62 <block dr="5." dz="0.07955" material="AlGaAs"/>
63 <block dr="5." dz="0.06949" material="GaAs"/>
64 </stack>
65 <item zero="0">
66 <stack name="middle">
67 <block name="x1" dr="5." dz="0.06371" material="AlGaAs"/>
68 <shelf name="oxide-layer">
69 <block dr="4." dz="0.01593" material="AlAs"/><block dr="1." dz="0.01593" material="AlOx"/>
70 </shelf>
71 <block name="x" dr="5." dz="0.00000" material="AlGaAs"/>
72 <block dr="5." dz="0.13649" material="GaAs"/>
73 <shelf name="QW">
74 <block name="active" role="gain" dr="4." dz="0.00500" material="InGaAs"/><block dr="1." dz="0.00500" material="InGaAs"/>
75 </shelf>
76 <zero/>
77 <block dr="5." dz="0.13649" material="GaAs"/>
78 </stack>
79 </item>
80 <stack name="bottom-dbr" repeat="29">
81 <block dr="5." dz="0.07955" material="AlGaAs"/>
82 <block dr="5." dz="0.06949" material="GaAs"/>
83 </stack>
84 <block dr="5." dz="0.07955" material="AlGaAs"/>
85 </stack>
86 </cylindrical>
87 </geometry>
88 <solvers>
89 <optical name="bessel" lib="modal" solver="BesselCyl">
90 <geometry ref="vcsel"/>
91 <expansion domain="finite" lam0="980." k-scale="0.1" k-method="laguerre"/>
92 <mode emission="top"/>
93 <pml dist="20." factor="1-0j" size="2.0"/>
94 <interface object="QW"/>
95 </optical>
96 </solvers>
97 </plask>''')
98 self.solver = self.manager.solvers.bessel
99 self.solver.initialize()
100 self.profile = StepProfile(self.solver.geometry)
101 self.solver.inGain = self.profile.outGain
102 self.solver.lam0 = 980.
103 self.solver.size = 30
104
105 def testComputationsFinite(self):
106 self.solver.domain = 'finite'
107 m = self.solver.find_mode(980.1)
108 self.assertEqual(m, 0)
109 self.assertEqual(len(self.solver.modes), 1)
110 # self.assertAlmostEqual(self.solver.modes[m].lam.real, 979.59, 2)
111 self.assertAlmostEqual(self.solver.modes[m].lam.real, 979.587, 3)
112 self.assertAlmostEqual(self.solver.modes[m].lam.imag, -0.02077, 3)
113
114 # Test integration of the Pointing vector
115 R = 27.
116 n = 1000
117 dr = 1e-6*R / n
118 rr = linspace(0., R, n+1)
119 msh = mesh.Rectangular2D(rr, [self.solver.geometry.bbox.top + 1e-6])
120 E = self.solver.outLightE(m, msh).array[:,0,:]
121 H = self.solver.outLightH(m, msh).array[:,0,:]
122 P = 0.5 * real(E[:,1]*conj(H[:,0]) - E[:,0]*conj(H[:,1]))
123 self.assertAlmostEqual(2e3*pi * sum(1e-6*rr * P) * dr / self.solver.modes[m].power, 1.0, 3)
124
125 def testComputationsInfinite(self):
126 self.solver.domain = 'infinite'
127 self.solver.kscale = 0.05
128 m = self.solver.find_mode(979.0)
129 self.assertEqual(m, 0)
130 self.assertEqual(len(self.solver.modes), 1)
131 self.assertAlmostEqual(self.solver.modes[m].lam.real, 979.663, 3)
132 self.assertAlmostEqual(self.solver.modes[m].lam.imag, -0.02077, 3)
133
134 # Test integration of the Pointing vector
135 R = 27.
136 n = 1000
137 dr = 1e-6*R / n
138 rr = linspace(0., R, n+1)
139 msh = mesh.Rectangular2D(rr, [self.solver.geometry.bbox.top + 1e-6])
140 E = self.solver.outLightE(m, msh).array[:,0,:]
141 H = self.solver.outLightH(m, msh).array[:,0,:]
142 P = 0.5 * real(E[:,1]*conj(H[:,0]) - E[:,0]*conj(H[:,1]))
143 self.assertAlmostEqual(2e3*pi * sum(1e-6*rr * P) * dr / self.solver.modes[m].power, 1.0, 2)
144
145 def _integrals_test(self, domain, prec, right=None, nr=101, kscale=None):
146 self.solver.domain = domain
147 if kscale is not None:
148 self.solver.kscale = kscale
149 bbox = self.solver.geometry.get_object_bboxes(self.manager.geo.middle)[0]
150 if right is None:
151 right = bbox.right
152 msh = mesh.Rectangular2D(mesh.Regular(0., right, nr), mesh.Regular(bbox.bottom, bbox.top, 501))
153 m = self.solver.find_mode(979.0)
154 dr = msh.axis0[1] - msh.axis0[0]
155 dz = msh.axis1[1] - msh.axis1[0]
156 integral_mesh = msh.elements.mesh
157 rr, _ = meshgrid(integral_mesh.axis0, integral_mesh.axis1, indexing='ij')
158
159 E = self.solver.outLightE(integral_mesh).array
160 E2 = sum(real(E*conj(E)), -1)
161 EE0 = 0.5 * 2*pi * sum((rr * E2).ravel()) * dr * dz
162 EE1 = self.solver.integrateEE(bbox.bottom, bbox.top)
163 self.assertAlmostEqual(EE0 / EE1, 1., prec)
164
165 H = self.solver.outLightH(integral_mesh).array
166 H2 = sum(real(H*conj(H)), -1)
167 HH0 = 0.5 * 2*pi * sum((rr * H2).ravel()) * dr * dz
168 HH1 = self.solver.integrateHH(bbox.bottom, bbox.top)
169 self.assertAlmostEqual(HH0 / HH1, 1., prec)
170
171 def testIntegralsFinite(self):
172 self._integrals_test('finite', 2, kscale=0.1)
173
174 def testIntegralsInfinite(self):
175 self._integrals_test('infinite', 2, 20., 301, kscale=0.05)
176
177 def plot_determinant(self):
178 lams = linspace(979., 982., 201)
179 dets = self.solver.get_determinant(lam=lams, m=1)
180 plot(lams, abs(dets))
181 yscale('log')
182
183 def plot_field(self):
184 self.solver.find_mode(979.0, m=1)
185 print_log('result', self.solver.modes[0])
186 box = self.solver.geometry.bbox
187 msh = mesh.Rectangular2D(mesh.Regular(0., box.right, 101),
188 mesh.Regular(box.bottom, box.top, 1001))
189 field = self.solver.outLightE(msh)
190 mag = max(abs(field.array.ravel()))
191 scale = linspace(-mag, mag, 255)
192 figure()
193 plot_geometry(self.solver.geometry, color='k', alpha=0.15)
194 plot_field(field, scale, comp='r', cmap='bwr')
195 window_title("Er")
196 colorbar(use_gridspec=True)
197 tight_layout(pad=0.1)
198 figure()
199 plot_geometry(self.solver.geometry, color='k', alpha=0.15)
200 plot_field(field, scale, comp='p', cmap='bwr')
201 colorbar(use_gridspec=True)
202 window_title("Ep")
203 tight_layout(pad=0.1)
204 figure()
205 plot_geometry(self.solver.geometry, color='k', alpha=0.15)
206 plot_field(field, scale, comp='z', cmap='bwr')
207 colorbar(use_gridspec=True)
208 window_title("Ez")
209 tight_layout(pad=0.1)
210
211 figure()
212 plot_geometry(self.solver.geometry, color='w', alpha=0.15)
213 light = self.solver.outLightMagnitude(msh)
214 plot_field(light)
215 colorbar(use_gridspec=True)
216 window_title("Mag")
217 tight_layout(pad=0.1)
218
219 z = self.solver.geometry.get_object_bboxes(self.manager.geo.QW)[0].center.z
220 arr = light.array
221 r = msh.axis0[int(unravel_index(argmax(arr), arr.shape)[0])]
222 rmsh = mesh.Rectangular2D(linspace(0, box.right, 2001), [z])
223 zmsh = mesh.Rectangular2D([r], linspace(box.bottom, box.top, 10001))
224 figure()
225 plot_profile(self.solver.outLightMagnitude(rmsh))
226 window_title(u"Horizontal (z = {:.1f} µm".format(z))
227 tight_layout(pad=0.1)
228 figure()
229 plot_profile(self.solver.outLightMagnitude(zmsh), swap_axes=True)
230 window_title(u"Vertical (r = {:.1f} µm".format(r))
231 tight_layout(pad=0.1)
232
233
234if __name__ == "__main__":
235 vcsel = VCSEL('plot_field')
236 vcsel.setUp()
237
238 try:
239 vcsel.solver.domain = sys.argv[1]
240 except IndexError:
241 pass
242
243 # vcsel.plot_determinant()
244 vcsel.plot_field()
245 show()