forked from Vlasko/MPhys
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFres_4.py
137 lines (101 loc) · 3.24 KB
/
Fres_4.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
from numpy import cos, inf, zeros, array, exp, conj, nan, isnan, pi, sin
import numpy as np
import scipy as sp
import time
start_time = time.time()
cos_theta0 = 1
# at first set theta1 to 0
cos_theta1 = 1
# when theta1 is 0, then so too is theta1
cos_theta2 = 1
# when theta1 is 0, then so too is theta2
n_0 = 1
# for air the refractive index is ~1
# refractive index of ITO glass
n_2 = 1.902
# code below can be used to choose between a substrate of Copper or Glass
# material = str(input("Is the substrate Copper or Glass? "))
#
# if material == "Copper":
# n_2 = 0.637
# elif material == "Glass":
# n_2 = 1.902
# else:
# print("invalid input")
# exit()
# choose between Cu or ITO glass refractive indecies
print(n_2)
# we want n to vary between 0 and 3 with a step of 0.1
steps = int(3/0.1) + 1
n = np.empty((steps))
filler = np.arange(0,3.1,0.1)
index_n = np.arange(n.size)
np.put(n,index_n,filler)
k = np.empty((steps))
index_k = np.arange(k.size)
np.put(k,index_k,filler)
import numpy as np
n = np.array(n)
k = np.array(k).reshape((-1, 1))
n_com = n + k.repeat(len(n), 1) * 1j
print (n_com)
print ( "time", time.time() - start_time)
r_01 = (cos_theta0 - n_com*cos_theta1)/(cos_theta0 + n_com*cos_theta1)
t_01 = (2*cos_theta0)/(cos_theta0 + n_com*cos_theta1)
# print(r_01)
# print(t_01)
r_12 = (n_com*cos_theta1 - n_2*cos_theta2)/(n_com*cos_theta1 + n_2*cos_theta2)
t_12 = (2*n_com*cos_theta1)/(n_com*cos_theta1 + n_2*cos_theta2)
# print(r_12)
# print(t_12)
lamda = 500 # in units nm
d = 200 # in units nm
beta = ((2*np.pi)/lamda)*n_com*d*cos_theta1
# print("Beta is equal to {}".format(beta))
z = 0 + 1j
block = np.exp(2*z*beta)
# not sure if this is defo doing what we hope
# but it is cycling through the set because we do have
# values at each "step"
# n is what we are changing
# should I be running over both i and j?
for i in r_01, r_12, block:
for j in i:
r_film = (r_01 + (r_12*block))/(1 + r_01*r_12*block)
# print r_film
print ( "time 2", time.time() - start_time)
# find r_film for all value of n contained in the array
for i in t_01, t_12, block:
for j in i:
t_film = (t_01*t_12*block)/(1 + r_01*r_12*block)
# find the absolute value of this and then square it
ABS_R = np.absolute(r_film)
# print(ABS_R)
R = np.power(ABS_R,2)
# find the value of R which is equal to square of absolute value of r
ABS_T = np.absolute(t_film)
# print(ABS_T)
T = n_2*np.power(ABS_T,2)
# find the value of T which is equal to square of absolute value of t
A = 1 - T - R
print ( "time 3", time.time() - start_time)
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(n,k,T, rstride=2, cstride=2)
plt.savefig('/Users/luka/Documents/University/MPhys/Theory/n&k_plot_T.jpg')
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(n,k,R, rstride=2, cstride=2)
ax.set_ylabel('k')
ax.set_xlabel('n')
plt.savefig('/Users/luka/Documents/University/MPhys/Theory/n&k_plot_R.jpg')
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(n,k,A, rstride=2, cstride=2)
ax.set_ylabel('k')
ax.set_xlabel('n')
plt.savefig('/Users/luka/Documents/University/MPhys/Theory/n&k_plot_A.jpg')