PLaSK library
Loading...
Searching...
No Matches
wire3d.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# coding=utf8
14
15import sys
16
17config.axes = 'zxy'
18config.log.level = 'detail'
19
20from optical.modal import Fourier3D
21
22@material.simple()
23class Glass(material.Material):
24 def Nr(self, wl, T, n): return 1.3
25
26solver = Fourier3D("solver")
27
28
29h_top = 0.1
30L = 6.0
31depth = 1.0
32
33
34symmetric = True
35
36interp = 'fourier'
37#interp = 'linear'
38
39
40provider = solver.outLightE
41#provider = solver.outLightH
42#provider = solver.outLightMagnitude
43comp = 1
44
45
46solver.size = 1, 32
47solver.smooth = 1e-3
48
49solver.vpml.dist = 1.00
50solver.vpml.size = 0.50
51solver.vpml.factor = 1-2j
52
53solver.pmls.tran.dist = L/2 - 0.25
54solver.pmls.tran.size = 0.5
55solver.pmls.tran.factor = 1-2j
56print_log('data', 'Tran. PML:', solver.pmls.tran)
57print_log('data', 'Vert. PML:', solver.vpml)
58
59
60solver.root.method = 'broyden'
61
62rect_top = geometry.Cuboid(depth, 0.25 if symmetric else 0.5, h_top, 'Glass')
63rect_bottom = geometry.Cuboid(depth, 0.25 if symmetric else 0.5, 1.5-h_top, 'Glass')
64
65stack = geometry.Stack3D(shift=-rect_bottom.height, back=0, **({'left':0} if symmetric else {'xcenter':0}))
66stack.prepend(rect_top)
67stack.prepend(rect_bottom)
68
69spacer = geometry.Cuboid(depth, L/2 if symmetric else L, 2., None)
70#stack.prepend(spacer)
71
72space = geometry.Cartesian3D(stack,
73 back='periodic', front='periodic',
74 left='mirror' if symmetric else 'air', right='air',
75 bottom='air', top='air')
76
77solver.geometry = space
78solver.set_interface(rect_top)
79solver.wavelength = 1000.
80
81def compute(sym):
82
83 if symmetric:
84 solver.symmetry = None, sym
85 else:
86 solver.symmetry = None, None
87
88 start = {'Ez': 1.1433, 'Ex': 1.115}[sym]
89
90 #nn = linspace(1.001, 1.299, 200)
91 #dets = solver.get_determinant(klong=nn*solver.k0)
92 #figure()
93 #plot(nn, abs(dets), 'r', label='PLaSK')
94 #xlabel('neff')
95 #legend(loc='best')
96 #gcf().canvas.set_window_title('Determinant [{}]'.format(sym))
97 #show()
98 #sys.exit(0)
99
100 mod = solver.find_mode(klong=start*solver.k0)
101
102 solver.modes[mod].power = 1.
103 print_log('result', f"neff: {solver.modes[mod].klong/solver.modes[mod].k0}", solver.modes[mod])
104
105 lx = (solver.pmls.tran.size + L/2.)
106 ly_b = - rect_bottom.height - solver.vpml.dist - solver.vpml.size
107 ly_t = rect_top.height + solver.vpml.dist + solver.vpml.size
108
109 XX = linspace(-lx, lx, 201)
110 YY = linspace(ly_b, ly_t, 201)
111 ZZ = [0.]
112
113 h = h_top - stack.bbox.height / 2
114
115 figure()
116 msh = mesh.Rectangular3D(ZZ, XX, YY)
117 field = solver.outLightMagnitude(mod, msh, interp)
118 #field = Data(sum(provider(mod, msh, interp).array[0,:,:,:]**2, 2), msh)
119 im = plot_field(field, plane='xy')
120 colorbar(im, use_gridspec=True)
121 plot_geometry(solver.geometry, 'w', mirror=True, plane='xy')
122 axvline(0, color='w', ls=':')
123 axhline(h, color='w', ls=':')
124 gcf().canvas.set_window_title('2D intensity [{}]'.format(sym))
125 tight_layout()
126
127 figure()
128 msh = mesh.Rectangular3D(ZZ, [0.], YY)
129 field = provider(mod, msh, interp).array[0,0,:,comp]
130 #plot(field.real, YY)
131 #plot(field.imag, YY)
132 plot(abs(field), YY)
133 bboxes = space.get_leafs_bboxes()
134 lines = set()
135 for box in bboxes:
136 lines.add(box.lower.y)
137 lines.add(box.upper.y)
138 for line in lines:
139 axhline(line, ls=':', color='k')
140 axhline(space.bbox.lower.y-solver.vpml.dist, ls=':', color='y')
141 axhline(space.bbox.upper.y+solver.vpml.dist, ls=':', color='y')
142 gcf().canvas.set_window_title("E{} (y) [{}]".format(config.axes[comp], sym))
143 ylim(min(YY), max(YY));
144 tight_layout()
145
146 figure()
147 msh = mesh.Rectangular3D(ZZ, XX, [h])
148 field = provider(mod, msh, interp).array[0,:,0,comp]
149 #plot(XX, field.real)
150 #plot(XX, field.imag)
151 plot(XX, abs(field))
152 xlim(XX[0],XX[-1])
153 bboxes = space.get_leafs_bboxes()
154 lines = set()
155 for box in bboxes:
156 lines.add(box.lower.x)
157 lines.add(box.upper.x)
158 for line in lines:
159 axvline(line, ls=':', color='k')
160 gcf().canvas.set_window_title("E{} (x) [{}]".format(config.axes[comp], sym))
161 tight_layout()
162
163 box = space.bbox
164 integral_mesh = mesh.Rectangular3D([-1, 1],
165 #mesh.Regular(-box.right if sym else box.left, box.right, 20001),
166 mesh.Regular(-lx, lx, 1001),
167 mesh.Regular(box.bottom, box.top, 1001))
168 dx, dy = integral_mesh.axis1[1] - integral_mesh.axis1[0], integral_mesh.axis2[1] - integral_mesh.axis2[0]
169 integral_mesh = integral_mesh.elements.mesh
170
171 E = solver.outLightE(mod, integral_mesh).array
172 E2 = sum(real(E*conj(E)), -1)
173 print_log('result', "E:", 0.5 * sum(E2) * dx * dy / solver.integrateEE(mod, box.bottom, box.top))
174
175 H = solver.outLightH(mod, integral_mesh).array
176 H2 = sum(real(H*conj(H)), -1)
177 print_log('result', "H:", 0.5 * sum(H2) * dx * dy / solver.integrateHH(mod, box.bottom, box.top))
178
179
180compute('Ez')
181compute('Ex')
182
183for mode in solver.modes:
184 print_log('important', f"neff: {mode.klong/mode.k0}", mode)
185
186show()