@@ -25,16 +25,6 @@ def BuildP(phil):
2525
2626### this function will compute the derivative of the P matrix for a given
2727### layer l with respect to the thickness of layer l
28- def Build_dP_ds (kzl , dl ):
29- P = np .zeros ((2 ,2 ),dtype = complex )
30- ci = 0 + 1j
31- a = - 1 * ci * kzl * dl
32- b = ci * kzl * dl
33- P [0 ][1 ] = 0 + 0j
34- P [1 ][0 ] = 0 + 0j
35- P [0 ][0 ] = - ci * kzl * np .exp (a )
36- P [1 ][1 ] = ci * kzl * np .exp (b )
37- return P
3828
3929def BuildD (nl , ctheta ,pol ):
4030 D = np .zeros ((2 ,2 ),dtype = complex )
@@ -172,151 +162,7 @@ def tmm(k0, theta0, pol, nA, tA):
172162 return M
173163
174164
175- ### analytically differentiates the transfer matrix
176- ### with respect to the thickness of a layer li to be specified by
177- ### the user!
178- def d_tmm_dsi (li , k0 , theta0 , pol , nA , tA ):
179- t1 = np .zeros ((2 ,2 ),dtype = complex )
180- t2 = np .zeros ((2 ,2 ),dtype = complex )
181- Dl = np .zeros ((2 ,2 ),dtype = complex )
182- Dli = np .zeros ((2 ,2 ),dtype = complex )
183- Pl = np .zeros ((2 ,2 ),dtype = complex )
184- M = np .zeros ((2 ,2 ),dtype = complex )
185-
186- D1 = BuildD (nA [0 ], np .cos (theta0 ), pol )
187- ### Note it is actually faster to invert the 2x2 matrix
188- ### "By Hand" than it is to use linalg.inv
189- ### and this inv step seems to be the bottleneck for the TMM function
190- tmp = D1 [0 ,0 ]* D1 [1 ,1 ]- D1 [0 ,1 ]* D1 [1 ,0 ]
191- det = 1 / tmp
192- M [0 ,0 ] = det * D1 [1 ,1 ]
193- M [0 ,1 ] = - det * D1 [0 ,1 ]
194- M [1 ,0 ] = - det * D1 [1 ,0 ]
195- M [1 ,1 ] = det * D1 [0 ,0 ]
196- #D1i = inv(D1)
197- #print("D1i is ")
198- #print(D1i)
199-
200-
201- ### This is the number of layers in the structure
202- L = len (nA )
203- kz = np .zeros (L ,dtype = complex )
204- phil = np .zeros (L ,dtype = complex )
205- ctheta = np .zeros (L ,dtype = complex )
206- theta = np .zeros (L ,dtype = complex )
207-
208- ### since kx is conserved through all layers, just compute it
209- ### in the upper layer (layer 1), for which you already known
210- ### the angle of incidence
211- kx = nA [0 ]* k0 * np .sin (theta0 )
212- kz [0 ] = np .sqrt ((nA [0 ]* k0 )** 2 - kx ** 2 )
213- kz [L - 1 ] = np .sqrt ((nA [L - 1 ]* k0 )** 2 - kx ** 2 )
214-
215- ### keeping consistent with K-R excitation
216- if np .real (kz [0 ])< 0 :
217- kz [0 ] = - 1 * kz [0 ]
218- if np .imag (kz [L - 1 ])< 0 :
219- kz [L - 1 ] = - 1 * kz [L - 1 ]
220-
221- ''' For dM/dsi we have three cases:
222- we compute contributions to M as normal for layers 2-li-1
223- we compute dPi/dsi for layer li
224- we compute contributions to M as normal for layers li+1 - L-1
225- '''
226-
227- ''' layers 2 - li-1 '''
228- for i in range (1 ,li ):
229- kz [i ] = np .sqrt ((nA [i ]* k0 )** 2 - kx ** 2 )
230- if np .imag (kz [i ])< 0 :
231- kz [i ] = - 1 * kz [i ]
232-
233- ctheta [i ] = kz [i ]/ (nA [i ]* k0 )
234- theta [i ] = np .arccos (ctheta [i ])
235-
236- phil [i ] = kz [i ]* tA [i ]
237-
238- Dl = BuildD (nA [i ],ctheta [i ], pol )
239- ## Invert Dl
240- tmp = Dl [0 ,0 ]* Dl [1 ,1 ]- Dl [0 ,1 ]* Dl [1 ,0 ]
241- det = 1 / tmp
242- Dli [0 ,0 ] = det * Dl [1 ,1 ]
243- Dli [0 ,1 ] = - det * Dl [0 ,1 ]
244- Dli [1 ,0 ] = - det * Dl [1 ,0 ]
245- Dli [1 ,1 ] = det * Dl [0 ,0 ]
246- #Dli = inv(Dl)
247- ## form Pl
248- Pl = BuildP (phil [i ])
249-
250- t1 = np .matmul (M ,Dl )
251- t2 = np .matmul (t1 ,Pl )
252- M = np .matmul (t2 ,Dli )
253-
254- ''' layer li '''
255- kz [li ] = np .sqrt ((nA [li ]* k0 )** 2 - kx ** 2 )
256- if np .imag (kz [li ])< 0 :
257- kz [li ] = - 1 * kz [li ]
258-
259- ctheta [li ] = kz [li ]/ (nA [li ]* k0 )
260- theta [li ] = np .arccos (ctheta [li ])
261-
262- Dl = BuildD (nA [li ], ctheta [li ], pol )
263- tmp = Dl [0 ,0 ]* Dl [1 ,1 ]- Dl [0 ,1 ]* Dl [1 ,0 ]
264- det = 1 / tmp
265- Dli [0 ,0 ] = det * Dl [1 ,1 ]
266- Dli [0 ,1 ] = - det * Dl [0 ,1 ]
267- Dli [1 ,0 ] = - det * Dl [1 ,0 ]
268- Dli [1 ,1 ] = det * Dl [0 ,0 ]
269-
270- Pl = Build_dP_ds (kz [li ], tA [li ])
271-
272- t1 = np .matmul (M , Dl )
273- t2 = np .matmul (t1 , Pl )
274- M = np .matmul (t2 , Dli )
275-
276-
277- ''' layers li+1 - L-1 '''
278- for i in range (li + 1 ,(L - 1 )):
279- kz [i ] = np .sqrt ((nA [i ]* k0 )** 2 - kx ** 2 )
280- if np .imag (kz [i ])< 0 :
281- kz [i ] = - 1 * kz [i ]
282-
283- ctheta [i ] = kz [i ]/ (nA [i ]* k0 )
284- theta [i ] = np .arccos (ctheta [i ])
285-
286- phil [i ] = kz [i ]* tA [i ]
287-
288- Dl = BuildD (nA [i ],ctheta [i ], pol )
289- ## Invert Dl
290- tmp = Dl [0 ,0 ]* Dl [1 ,1 ]- Dl [0 ,1 ]* Dl [1 ,0 ]
291- det = 1 / tmp
292- Dli [0 ,0 ] = det * Dl [1 ,1 ]
293- Dli [0 ,1 ] = - det * Dl [0 ,1 ]
294- Dli [1 ,0 ] = - det * Dl [1 ,0 ]
295- Dli [1 ,1 ] = det * Dl [0 ,0 ]
296- #Dli = inv(Dl)
297- ## form Pl
298- Pl = BuildP (phil [i ])
299-
300- t1 = np .matmul (M ,Dl )
301- t2 = np .matmul (t1 ,Pl )
302- M = np .matmul (t2 ,Dli )
303-
304- ### M is now the product of D_1^-1 .... D_l-1^-1... just need to
305- ### compute D_L and multiply M*D_L
306- kz [L - 1 ] = np .sqrt ((nA [L - 1 ]* k0 )** 2 - kx ** 2 )
307- ctheta [L - 1 ]= kz [L - 1 ]/ (nA [L - 1 ]* k0 )
308- DL = BuildD (nA [L - 1 ], ctheta [L - 1 ], pol )
309- t1 = np .matmul (M ,DL )
310- ### going to create a dictionary called M which will
311- ### contain the matrix elements of M as well as
312- ### other important quantities like incoming and outgoing angles
313- dM_ds = {"dM11_ds" : t1 [0 ,0 ],
314- "dM12_ds" : t1 [0 ,1 ],
315- "dM21_ds" : t1 [1 ,0 ],
316- "dM22_ds" : t1 [1 ,1 ]
317- }
318165
319- return dM_ds
320166### This function will take a thickness array
321167### and a scalar thickness and will return
322168### the index of the layer in which you are in
0 commit comments