33import numpy
44import numpy .linalg
55
6+
67def build_anisotropic_permittivity (e1 , e2 , e3 ):
78 eps = numpy .zeros ([3 , 3 ], dtype = numpy .complex128 )
89 eps [0 , 0 ] = e1
910 eps [1 , 1 ] = e2
1011 eps [2 , 2 ] = e3
1112 return eps
1213
14+
1315def build_rotation_matrix (theta , phi , psi ):
14- # Rx = numpy.asarray(
15- # [1, 0, 0],
16- # [0, numpy.cos(theta), -numpy.sin(theta)],
17- # [0, numpy.sin(theta), numpy.cos(theta)]
18- # )
19- #
20- # Ry = numpy.asarray(
21- # [numpy.cos(phi), 0, numpy.sin(phi)],
22- # [0, 1, 0],
23- # [-numpy.sin(phi), numpy.cos(phi)]
24- # )
25- #
26- # Rz = numpy.asarray(
27- # [numpy.cos(psi), -numpy.sin(psi), 0],
28- # [numpy.sin(psi), numpy.cos(psi), 0],
29- # [0, 0, 1]
30- # )
31- #
32- # R = numpy.dot(Rz, numpy.dot(Ry, Rx))
3316 R = numpy .asarray ([
3417 [numpy .cos (psi )* numpy .cos (phi ) - numpy .cos (theta )* numpy .sin (phi )* numpy .sin (psi ),
3518 - numpy .sin (psi )* numpy .cos (phi ) - numpy .cos (theta )* numpy .sin (phi )* numpy .cos (psi ),
@@ -42,6 +25,7 @@ def build_rotation_matrix(theta, phi, psi):
4225
4326 return R
4427
28+
4529def build_general_permittivity_matrix (e1 , e2 , e3 , theta , phi , psi ):
4630 E = build_anisotropic_permittivity (e1 , e2 , e3 )
4731 R = build_rotation_matrix (theta , phi , psi )
@@ -54,131 +38,8 @@ def build_general_permittivity_matrix(e1, e2, e3, theta, phi, psi):
5438
5539 return Eps
5640
57- def build_polarization_vector (w , Eps , kx , ky , g , mu ):
58-
59- e_zz = Eps [2 , 2 ]
60- e_yy = Eps [1 , 1 ]
61- e_xx = Eps [0 , 0 ]
62- e_xy = Eps [0 , 1 ]
63- e_yz = Eps [1 , 2 ]
64- e_xz = Eps [0 , 2 ]
65-
66- # TODO: Guaranteed bullshit?
67- p = numpy .asarray (
68- [(w ** 2 * e_yy - kx ** 2 - g ** 2 )* (w ** 2 * e_zz - kx ** 2 - ky ** 2 ) - (w ** 2 * e_yz + ky * g )** 2 ,
69- (w ** 2 * e_yz + ky * g )* (w ** 2 * e_xz + kx * g ) - (w ** 2 * e_xy + kx * ky )* (w ** 2 * e_zz - kx ** 2 - ky ** 2 ),
70- (w ** 2 * e_xy + kx * ky )* (w ** 2 * e_yz + ky * g ) - (w ** 2 * e_xz + kx * g )* (w ** 2 * e_yy - kx ** 2 - g ** 2 )],
71- dtype = numpy .complex128 )
72- p = numpy .divide (p , numpy .sqrt (numpy .dot (p , p )))
73- return p
74-
75- def build_anisotropic_layer_matrix (e1 , e2 , e3 , theta , phi , psi , w , kx , ky , d ):
76- #
77- # Build epsilon matrix
78- #
79- Eps = build_general_permittivity_matrix (e1 , e2 , e3 , theta , phi , psi )
80-
81- e_zz = Eps [2 , 2 ]
82- e_yy = Eps [1 , 1 ]
83- e_xx = Eps [0 , 0 ]
84- e_xy = Eps [0 , 1 ]
85- e_yz = Eps [1 , 2 ]
86- e_xz = Eps [0 , 2 ]
87- # should be symmetric, so we don't need e_zy, e_zx, e_yx
88-
89- #
90- # Solve for eigenmodes
91- # Solve quartic equation A*kz**4 + B*kz**3 + C*kz**2 + D*kz + E
92- # Coefficients from Determinant of matrix
93- # [(w/c)**2 * E + inner(k, k) * I + outer(k, k)]
94- a = w ** 2 * e_zz # *kz**4
95-
96- b = w ** 2 * (2 * e_xz * kx + 2 * e_yz * ky ) # *kz**3
97-
98- c = w ** 2 * ( - e_xx * e_zz * w ** 2 + e_xx * kx ** 2 + 2 * e_xy * kx * ky + e_xz ** 2 * w ** 2 - e_yy * e_zz * w ** 2 + e_yy * ky ** 2
99- + e_yz ** 2 * w ** 2 + e_zz * kx ** 2 + e_zz * ky ** 2
100- ) # *kz**2
101-
102- d = w ** 2 * ( - 2 * e_xx * e_yz * w ** 2 * ky + 2 * e_xy * e_xz * w ** 2 * ky + 2 * e_xy * e_yz ** w ** 2 * kx - 2 * e_xz * e_yy * w ** 2 * kx + 2 * e_xz * kx ** 3
103- + 2 * e_xz * kx * ky ** 2 + 2 * e_yz * kx ** 2 * ky + 2 * e_yz * ky ** 3
104- ) # *kz**1
105-
106- e = w ** 2 * (e_xx * e_yy * e_zz * w ** 4 - e_xx * e_yy * w ** 2 * kx ** 2 - e_xx * e_yy * w ** 2 * ky ** 2 - e_xx * e_yz ** 2 * w ** 4
107- - e_xx * e_zz * w ** 2 * kx ** 2 + e_xx * kx ** 4 + e_xx * kx ** 2 * ky ** 2 - e_xy ** 2 * e_zz * w ** 4 + e_xy ** 2 * w ** 2 * kx ** 2
108- + e_xy ** 2 * w ** 2 * ky ** 2 + 2 * e_xy * e_xz * e_yz * w ** 4
109- - 2 * e_xy * e_zz * w ** 2 * kx * ky + 2 * e_xy * kx ** 3 * ky + 2 * e_xy * kx * ky ** 3 - e_xz ** 2 * e_yy * w ** 4
110- + e_xz ** 2 * w ** 2 * kx ** 2 + 2 * e_xz * e_yz * w ** 2 * kx * ky - e_yy * e_zz * w ** 2 * ky ** 2 + e_yy * kx ** 2 * ky ** 2
111- + e_yy * ky ** 4 + e_yz ** 2 * w ** 2 * ky ** 2
112- ) # *kz**0
113- coeffs = [a , b , c , d , e ]
114- gamma = numpy .roots (coeffs )
115-
116- #TODO: Hack needs fixing
117- tmp = gamma [2 ]
118- gamma [2 ] = gamma [1 ]
119- gamma [1 ] = tmp
120-
121- #
122- # Build polarization vectors
123- #
124- mu = 1
125- c = 1 # m/c
126-
127- p = [build_polarization_vector (w , Eps , kx , ky , g , mu ) for g in gamma ]
128- print ("P:" , numpy .real (p ))
129- q = [(c / (w * mu ))* numpy .cross ([kx , ky , gi ], pi ) for gi , pi in zip (gamma , p )]
130- print ("Q:" , numpy .real (q ))
131-
132- #return p, gamma, Eps
133-
134- #
135- # Tests
136- #
137- #vec = numpy.dot(Eps, p[0])
138- #print( numpy.dot(vec, [kx, ky, gamma[0]]) )
139- #
140- #
141- #
142-
143- #
144- # Build boundary transition matrix
145- #
146- D = numpy .asarray (
147- [
148- [p [0 ][0 ], p [1 ][0 ], p [2 ][0 ], p [3 ][0 ]],
149- [q [0 ][1 ], q [1 ][1 ], q [2 ][1 ], q [3 ][1 ]],
150- [p [0 ][1 ], p [1 ][1 ], p [2 ][1 ], p [3 ][1 ]],
151- [q [0 ][0 ], q [1 ][0 ], q [2 ][0 ], q [3 ][0 ]]
152- ], dtype = numpy .complex128
153- )
154- print ("D:" , numpy .real (D ))
155-
156- #
157- # Build propagation matrix
158- #
159- P = numpy .asarray (
160- [
161- [numpy .exp (- 1j * gamma [0 ]* d ), 0 , 0 , 0 ],
162- [0 , numpy .exp (- 1j * gamma [1 ]* d ), 0 , 0 ],
163- [0 , 0 , numpy .exp (- 1j * gamma [2 ]* d ), 0 ],
164- [0 , 0 , 0 , numpy .exp (- 1j * gamma [3 ]* d )]
165- ], dtype = numpy .complex128
166- )
167-
168- #
169- # Multiply matricies
170- #
171- LayerMatrix = numpy .dot (D , numpy .dot (P , numpy .linalg .inv (D )))
172- #return LayerMatrix
173- return D
174-
17541
17642def solve_transfer_matrix (M ):
177-
178- denom = M [0 , 0 ]* M [2 , 2 ] - M [0 , 2 ]* M [2 , 0 ]
179- #print(denom)
180-
181-
18243 #
18344 # Components 3,4 - S
18445 # Components 1,2 - P Was mixed up until 11.08.2016 (switched s and p)
0 commit comments