PLaSK library
Loading...
Searching...
No Matches
fresnel.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 sys
15import unittest
16
17from numpy import *
18
19from plask import *
20from plask import material, geometry, mesh
21from optical.modal import Fourier2D, Fourier3D
22
23plask.config.axes = 'xyz'
24
25angles = linspace(0.0, 90.0, 11)
26
27nr = material.GaAs().nr(1000)
28print("nr = {}".format(nr))
29
30class Refl:
31
32 def __init__(self, solver, direction, polarization):
33 self.solver = solver
34 self.solver.lam0 = 1000.
35 self.direction = direction
36 self.polarization = polarization
37 solver.klong = solver.ktran = 0.
38
39 def __call__(self, a):
40 setattr(self.solver, 'k'+self.direction, 2*pi * sin(a*pi/180.))
41 return self.solver.compute_reflectivity(1000., 'top', self.polarization), \
42 self.solver.compute_transmittivity(1000., 'top', self.polarization)
43
44 def __enter__(self):
45 return self
46
47 def __exit__(self, *args, **kwargs):
48 pass
49
50
51def R_TE(angle):
52 angle = angle * pi/180
53 cos1 = cos(angle)
54 cos2 = sqrt(1. - (sin(angle)/nr)**2)
55 return 100. * ((cos1 - nr*cos2) / (cos1 + nr*cos2))**2
56
57def R_TM(angle):
58 angle = angle * pi/180
59 cos1 = cos(angle)
60 cos2 = sqrt(1. - (sin(angle)/nr)**2)
61 return 100. * ((cos2 - nr*cos1) / (cos2 + nr*cos1))**2
62
63def T_TE(angle):
64 return 100. - R_TE(angle)
65
66def T_TM(angle):
67 return 100. - R_TM(angle)
68
69
70def show_plots(solver, direction, title=None, sep=False):
71 solver.klong = solver.ktran = 0.
72 if sep: solver.polarization = 'El'
73 with Refl(solver, direction, 'El') as refl:
74 Rlong = array([refl(a) for a in angles]).T
75 if sep: solver.polarization = 'Et'
76 with Refl(solver, direction, 'Et') as refl:
77 Rtran = array([refl(a) for a in angles]).T
78 f = figure()
79 plot(angles, Rlong[0], label='R long')
80 plot(angles, Rtran[0], label='R tran')
81 plot(angles, Rlong[1], label='T long')
82 plot(angles, Rtran[1], label='T tran')
83 xlabel('angle [deg]')
84 ylabel('reflectivity/transmittivity [%]')
85 ylim(0,100)
86 legend(loc='best')
87 plask.window_title((direction if title is None else title).title())
88
89
90class FresnelTest(unittest.TestCase):
91
92 def testFresnel2D(self):
93 top = geometry.Rectangle(1, 0.5, 'air')
94 bottom = geometry.Rectangle(1, 0.5, 'GaAs')
95 stack = geometry.Stack2D()
96 stack.prepend(top)
97 stack.prepend(bottom)
98 geo = geometry.Cartesian2D(stack, left='periodic', right='periodic', top='extend', bottom='extend')
99
100 solver = Fourier2D()
101 solver.geometry = geo
102 solver.size = 3
103
104 solver.polarization = 'El'
105 with Refl(solver, 'tran', 'El') as refl:
106 for a in angles:
107 R, T = refl(a)
108 self.assertAlmostEqual(R, R_TE(a), 3)
109 self.assertAlmostEqual(T, T_TE(a), 3)
110
111 solver.polarization = 'Et'
112 with Refl(solver, 'tran', 'Et') as refl:
113 for a in angles:
114 R, T = refl(a)
115 self.assertAlmostEqual(R, R_TM(a), 3)
116 self.assertAlmostEqual(T, T_TM(a), 3)
117
118 if __name__ == '__main__':
119 show_plots(solver, 'tran', 'Separated Tran 2D', True)
120
121 solver.polarization = None
122 #solver.initialize()
123 with Refl(solver, 'long', 'TE') as refl:
124 for a in angles:
125 R, T = refl(a)
126 self.assertAlmostEqual(R, R_TM(a), 3)
127 self.assertAlmostEqual(T, T_TM(a), 3)
128
129 with Refl(solver, 'long', 'TM') as refl:
130 for a in angles:
131 R, T = refl(a)
132 self.assertAlmostEqual(R, R_TE(a), 3)
133 self.assertAlmostEqual(T, T_TE(a), 3)
134
135 with Refl(solver, 'tran', 'TE') as refl:
136 for a in angles:
137 R, T = refl(a)
138 self.assertAlmostEqual(R, R_TE(a), 3)
139 self.assertAlmostEqual(T, T_TE(a), 3)
140
141 with Refl(solver, 'tran', 'TM') as refl:
142 for a in angles:
143 R, T = refl(a)
144 self.assertAlmostEqual(R, R_TM(a), 3)
145 self.assertAlmostEqual(T, T_TM(a), 3)
146
147 if __name__ == '__main__':
148 show_plots(solver, 'long', 'Long 2D')
149 show_plots(solver, 'tran', 'Tran 2D')
150
151
152 def testFresnel3D(self):
153 top = geometry.Cuboid(1, 1, 0.5, 'air')
154 bottom = geometry.Cuboid(1, 1, 0.5, 'GaAs')
155 stack = geometry.Stack3D()
156 stack.prepend(top)
157 stack.prepend(bottom)
158 geo = geometry.Cartesian3D(stack, left='periodic', right='periodic', back='periodic', front='periodic', top='extend', bottom='extend')
159
160 solver = Fourier3D()
161 solver.geometry = geo
162 solver.size = 1,1
163
164 with Refl(solver, 'long', 'El') as refl:
165 for a in angles:
166 R, T = refl(a)
167 self.assertAlmostEqual(R, R_TM(a), 3)
168 self.assertAlmostEqual(T, T_TM(a), 3)
169
170 with Refl(solver, 'long', 'Et') as refl:
171 for a in angles:
172 R, T = refl(a)
173 self.assertAlmostEqual(R, R_TE(a), 3)
174 self.assertAlmostEqual(T, T_TE(a), 3)
175
176 with Refl(solver, 'tran', 'El') as refl:
177 for a in angles:
178 R, T = refl(a)
179 self.assertAlmostEqual(R, R_TE(a), 3)
180 self.assertAlmostEqual(T, T_TE(a), 3)
181
182 with Refl(solver, 'tran', 'Et') as refl:
183 for a in angles:
184 R, T = refl(a)
185 self.assertAlmostEqual(R, R_TM(a), 3)
186 self.assertAlmostEqual(T, T_TM(a), 3)
187
188 if __name__ == '__main__':
189 show_plots(solver, 'long', 'Long 3D')
190 show_plots(solver, 'tran', 'Tran 3D')
191
192
193if __name__ == '__main__':
194 test = unittest.main(exit=False)
195 show()
196 sys.exit(not test.result.wasSuccessful())