-
Notifications
You must be signed in to change notification settings - Fork 93
/
examples.py
257 lines (238 loc) · 9.84 KB
/
examples.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
# -*- coding: utf-8 -*-
"""
Examples of plots and calculations using the tmm package.
"""
from __future__ import division, print_function, absolute_import
from .tmm_core import (coh_tmm, unpolarized_RT, ellips,
position_resolved, find_in_structure_with_inf)
from numpy import pi, linspace, inf, array
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
try:
import colorpy.illuminants
import colorpy.colormodels
from . import color
colors_were_imported = True
except ImportError:
# without colorpy, you can't run sample5(), but everything else is fine.
colors_were_imported = False
# "5 * degree" is 5 degrees expressed in radians
# "1.2 / degree" is 1.2 radians expressed in degrees
degree = pi/180
def sample1():
"""
Here's a thin non-absorbing layer, on top of a thick absorbing layer, with
air on both sides. Plotting reflected intensity versus wavenumber, at two
different incident angles.
"""
# list of layer thicknesses in nm
d_list = [inf, 100, 300, inf]
# list of refractive indices
n_list = [1, 2.2, 3.3+0.3j, 1]
# list of wavenumbers to plot in nm^-1
ks = linspace(0.0001, .01, num=400)
# initialize lists of y-values to plot
Rnorm = []
R45 = []
for k in ks:
# For normal incidence, s and p polarizations are identical.
# I arbitrarily decided to use 's'.
Rnorm.append(coh_tmm('s', n_list, d_list, 0, 1/k)['R'])
R45.append(unpolarized_RT(n_list, d_list, 45*degree, 1/k)['R'])
kcm = ks * 1e7 #ks in cm^-1 rather than nm^-1
plt.figure()
plt.plot(kcm, Rnorm, 'blue', kcm, R45, 'purple')
plt.xlabel('k (cm$^{-1}$)')
plt.ylabel('Fraction reflected')
plt.title('Reflection of unpolarized light at 0$^\circ$ incidence (blue), '
'45$^\circ$ (purple)')
plt.show()
def sample2():
"""
Here's the transmitted intensity versus wavelength through a single-layer
film which has some complicated wavelength-dependent index of refraction.
(I made these numbers up, but in real life they could be read out of a
graph / table published in the literature.) Air is on both sides of the
film, and the light is normally incident.
"""
#index of refraction of my material: wavelength in nm versus index.
material_nk_data = array([[200, 2.1+0.1j],
[300, 2.4+0.3j],
[400, 2.3+0.4j],
[500, 2.2+0.4j],
[750, 2.2+0.5j]])
material_nk_fn = interp1d(material_nk_data[:,0].real,
material_nk_data[:,1], kind='quadratic')
d_list = [inf, 300, inf] #in nm
lambda_list = linspace(200, 750, 400) #in nm
T_list = []
for lambda_vac in lambda_list:
n_list = [1, material_nk_fn(lambda_vac), 1]
T_list.append(coh_tmm('s', n_list, d_list, 0, lambda_vac)['T'])
plt.figure()
plt.plot(lambda_list, T_list)
plt.xlabel('Wavelength (nm)')
plt.ylabel('Fraction of power transmitted')
plt.title('Transmission at normal incidence')
plt.show()
def sample3():
"""
Here is a calculation of the psi and Delta parameters measured in
ellipsometry. This reproduces Fig. 1.14 in Handbook of Ellipsometry by
Tompkins, 2005.
"""
n_list = [1, 1.46, 3.87+0.02j]
ds = linspace(0, 1000, num=100) #in nm
psis = []
Deltas = []
for d in ds:
e_data = ellips(n_list, [inf, d, inf], 70*degree, 633) #in nm
psis.append(e_data['psi']/degree) # angle in degrees
Deltas.append(e_data['Delta']/degree) # angle in degrees
plt.figure()
plt.plot(ds, psis, ds, Deltas)
plt.xlabel('SiO2 thickness (nm)')
plt.ylabel('Ellipsometric angles (degrees)')
plt.title('Ellipsometric parameters for air/SiO2/Si, varying '
'SiO2 thickness.\n'
'@ 70$^\circ$, 633nm. '
'Should agree with Handbook of Ellipsometry Fig. 1.14')
plt.show()
def sample4():
"""
Here is an example where we plot absorption and Poynting vector
as a function of depth.
"""
d_list = [inf, 100, 300, inf] #in nm
n_list = [1, 2.2+0.2j, 3.3+0.3j, 1]
th_0 = pi/4
lam_vac = 400
pol = 'p'
coh_tmm_data = coh_tmm(pol, n_list, d_list, th_0, lam_vac)
ds = linspace(-50, 400, num=1000) #position in structure
poyn = []
absor = []
for d in ds:
layer, d_in_layer = find_in_structure_with_inf(d_list, d)
data = position_resolved(layer, d_in_layer, coh_tmm_data)
poyn.append(data['poyn'])
absor.append(data['absor'])
# convert data to numpy arrays for easy scaling in the plot
poyn = array(poyn)
absor = array(absor)
plt.figure()
plt.plot(ds, poyn, 'blue', ds, 200*absor, 'purple')
plt.xlabel('depth (nm)')
plt.ylabel('AU')
plt.title('Local absorption (purple), Poynting vector (blue)')
plt.show()
def sample5():
"""
Color calculations: What color is a air / thin SiO2 / Si wafer?
"""
if not colors_were_imported:
print('Colorpy was not detected (or perhaps an error occurred when',
'loading it). You cannot do color calculations, sorry!',
'Original version is at http://pypi.python.org/pypi/colorpy',
'A Python 3 compatible edit is at https://github.com/fish2000/ColorPy/')
return
# Crystalline silicon refractive index. Data from Palik via
# http://refractiveindex.info, I haven't checked it, but this is just for
# demonstration purposes anyway.
Si_n_data = [[400, 5.57 + 0.387j],
[450, 4.67 + 0.145j],
[500, 4.30 + 7.28e-2j],
[550, 4.08 + 4.06e-2j],
[600, 3.95 + 2.57e-2j],
[650, 3.85 + 1.64e-2j],
[700, 3.78 + 1.26e-2j]]
Si_n_data = array(Si_n_data)
Si_n_fn = interp1d(Si_n_data[:,0], Si_n_data[:,1], kind='linear')
# SiO2 refractive index (approximate): 1.46 regardless of wavelength
SiO2_n_fn = lambda wavelength : 1.46
# air refractive index
air_n_fn = lambda wavelength : 1
n_fn_list = [air_n_fn, SiO2_n_fn, Si_n_fn]
th_0 = 0
# Print the colors, and show plots, for the special case of 300nm-thick SiO2
d_list = [inf, 300, inf]
reflectances = color.calc_reflectances(n_fn_list, d_list, th_0)
illuminant = colorpy.illuminants.get_illuminant_D65()
spectrum = color.calc_spectrum(reflectances, illuminant)
color_dict = color.calc_color(spectrum)
print('air / 300nm SiO2 / Si --- rgb =', color_dict['rgb'], ', xyY =', color_dict['xyY'])
plt.figure()
color.plot_reflectances(reflectances,
title='air / 300nm SiO2 / Si -- '
'Fraction reflected at each wavelength')
plt.figure()
color.plot_spectrum(spectrum,
title='air / 300nm SiO2 / Si -- '
'Reflected spectrum under D65 illumination')
# Calculate irgb color (i.e. gamma-corrected sRGB display color rounded to
# integers 0-255) versus thickness of SiO2
max_SiO2_thickness = 600
SiO2_thickness_list = linspace(0, max_SiO2_thickness, num=80)
irgb_list = []
for SiO2_d in SiO2_thickness_list:
d_list = [inf, SiO2_d, inf]
reflectances = color.calc_reflectances(n_fn_list, d_list, th_0)
illuminant = colorpy.illuminants.get_illuminant_D65()
spectrum = color.calc_spectrum(reflectances, illuminant)
color_dict = color.calc_color(spectrum)
irgb_list.append(color_dict['irgb'])
# Plot those colors
print('Making color vs SiO2 thickness graph. Compare to (for example)')
print('http://www.htelabs.com/appnotes/sio2_color_chart_thermal_silicon_dioxide.htm')
plt.figure()
plt.plot([0, max_SiO2_thickness], [1, 1])
plt.xlim(0, max_SiO2_thickness)
plt.ylim(0, 1)
plt.xlabel('SiO2 thickness (nm)')
plt.yticks([])
plt.title('Air / SiO2 / Si color vs SiO2 thickness')
for i in range(len(SiO2_thickness_list)):
# One strip of each color, centered at x=SiO2_thickness_list[i]
if i == 0:
x0 = 0
else:
x0 = (SiO2_thickness_list[i] + SiO2_thickness_list[i-1]) / 2
if i == len(SiO2_thickness_list) - 1:
x1 = max_SiO2_thickness
else:
x1 = (SiO2_thickness_list[i] + SiO2_thickness_list[i+1]) / 2
y0 = 0
y1 = 1
poly_x = [x0, x1, x1, x0]
poly_y = [y0, y0, y1, y1]
color_string = colorpy.colormodels.irgb_string_from_irgb(irgb_list[i])
plt.fill(poly_x, poly_y, color_string, edgecolor=color_string)
plt.show()
def sample6():
"""
An example reflection plot with a surface plasmon resonance (SPR) dip.
Compare with http://doi.org/10.2320/matertrans.M2010003 ("Spectral and
Angular Responses of Surface Plasmon Resonance Based on the Kretschmann
Prism Configuration") Fig 6a
"""
# list of layer thicknesses in nm
d_list = [inf, 5, 30, inf]
# list of refractive indices
n_list = [1.517, 3.719+4.362j, 0.130+3.162j, 1]
# wavelength in nm
lam_vac = 633
# list of angles to plot
theta_list = linspace(30*degree, 60*degree, num=300)
# initialize lists of y-values to plot
Rp = []
for theta in theta_list:
Rp.append(coh_tmm('p', n_list, d_list, theta, lam_vac)['R'])
plt.figure()
plt.plot(theta_list/degree, Rp, 'blue')
plt.xlabel('theta (degree)')
plt.ylabel('Fraction reflected')
plt.xlim(30, 60)
plt.ylim(0, 1)
plt.title('Reflection of p-polarized light with Surface Plasmon Resonance\n'
'Compare with http://doi.org/10.2320/matertrans.M2010003 Fig 6a')
plt.show()