PLaSK library

Loading...
Searching...
No Matches
fields3d.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
14 import unittest
15
16 from numpy import *
17
18 from plask import *
19 from plask import material, geometry, mesh
20 from optical import modal
21
22 config.axes = 'zxy'
23
24@material.simple()
25 class Glass(material.Material):
26 def Nr(self, wl, T=300., n=0.): return 1.3
27
28
29solver = modal.Fourier3D("solver")
30 solver.ft = 'analytic'
31wire_stack = geometry.Stack2D()
32 wire_stack.append(geometry.Rectangle(0.75, 0.125, Glass()))
33rect = geometry.Rectangle(0.75, 0.375, Glass())
34 rect.role = 'interface'
35 wire_stack.append(rect)
36
37shelf = geometry.Shelf2D()
38 shelf.append(wire_stack)
39 shelf.append(geometry.Rectangle(0.15, 0.50, 'air'))
40
41extrusion = geometry.Extrusion(shelf, 1.)
42 solver.geometry = geometry.Cartesian3D(extrusion,
43 left="mirror", right="periodic",
44 back='mirror', front='periodic')
45
46 solver.dct = 1
47
48 solver.size = 0, 12
49
50start = 1.1216
51
52det = False
53
54 def plotFields(symmetry):
55 solver.lam = 1000.
56 solver.symmetry = None, symmetry
57
58 k0 = 2e3*pi / solver.lam
59
60 if det:
61 figure()
62 neffs = linspace(1.0, 1.2, 1001)
63 dets = solver.get_determinant(klong=k0*neffs)
64 plot(neffs, abs(dets))
65 window_title(f"det {symmetry}")
66
67 mn = solver.find_mode(klong=k0*start)
68 print_log('result', f"neff = {(solver.modes[mn].klong/k0).real}")
69
70 msh = mesh.Rectangular3D([0], mesh.Regular(-2., 2., 401), mesh.Regular(-0.4, 1.0, 1501))
71 figure(figsize=(14, 6))
72
73 E = solver.outLightE(mn, msh)
74 m = max(abs(E))
75 E /= m
76
77 ax_mag = subplot2grid((3,3), (0,1), rowspan=3)
78 plot_field(Data(sum(real(E*conj(E)),-1), E.mesh), plane='xy')
79 plot_geometry(solver.geometry, color='w', mirror=True, periods=2, lw=0.5, plane='xy')
80 xticks(arange(-2., 2.01, 0.5))
81 c = 2
82
83 subplot2grid((3,3), (0,2), rowspan=3, sharey=ax_mag)
84 msh1 = mesh.Rectangular3D([0], [0], msh.axis2)
85 E1 = solver.outLightE(mn, msh1)
86 E1 /= m
87 mr = max(abs(E.array[:,:,c].real).ravel())
88 mi = max(abs(E.array[:,:,c].imag).ravel())
89 if mi > mr: E1 = E1.imag
90 else: E1 = E1.real
91 axhline(0.0, ls='--', color='k', lw=0.5)
92 axhline(0.5, ls='--', color='k', lw=0.5)
93 plot_profile(E1, swap_axes=True, comp=c)
94 Ec = E1.array[0,0,:,c].copy()
95 y = array(msh.axis2)
96 Ec[(0.0 <= y) & (y <= 0.5)] *= 1.3**2
97 plot(Ec, y)
98 ylim(-0.4, 1.0)
99 window_title(f"Field {symmetry}")
100
101 levels = linspace(-1, 1, 16)
102 for c in range(3):
103 subplot2grid((3,3), (c, 0), sharey=ax_mag)
104 mr = max(abs(E.array[:,:,:,c].real).ravel())
105 mi = max(abs(E.array[:,:,:,c].imag).ravel())
106 print_log('data', f"Max E{config.axes[c]}: Re={mr} Im={mi}")
107 if mi > mr: Ec = E.imag
108 else: Ec = E.real
109 plot_field(Ec, levels, comp=c, cmap='bwr', plane='xy')
110 plot_geometry(solver.geometry, color='k', mirror=True, periods=2, lw=0.5, plane='xy')
111 xlim(msh.axis0[0], msh.axis0[-1])
112 xticks(arange(-2., 2.01, 0.5))
113
114 tight_layout()
115
116
117 plotFields('Htran')
118 plotFields(None)
119 show()