1 |
;;; "fresneleq.scm" Solve EM waves in parallel layers of dielectrics and metals -*-scheme-*- |
2 |
;;; Copyright (C) 2003, 2004 Aubrey Jaffer |
3 |
|
4 |
;;; This program is free software; you can redistribute it and/or modify |
5 |
;;; it under the terms of the GNU General Public License as published by |
6 |
;;; the Free Software Foundation, either version 3 of the License, or |
7 |
;;; (at your option) any later version. |
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 |
;;; |
14 |
;;; You should have received a copy of the GNU General Public License |
15 |
;;; along with this program; if not, write to the Free Software |
16 |
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA. |
17 |
|
18 |
;;; I can't find a free implementation of thin-film optical |
19 |
;;; calculations. So I rolled my own. |
20 |
|
21 |
;; http://people.csail.mit.edu/jaffer/FreeSnell |
22 |
|
23 |
;; Optics is all wavelength based; so this is also. These routines |
24 |
;; work out the complex voltages in the forward and reverse directions |
25 |
;; to find the transmitted and reflected amplitudes. This works only |
26 |
;; for intensities where the layers act linearly (superposition). |
27 |
;; Each layer 0:n has an index of refraction and height. |
28 |
|
29 |
;; Square the absolute value of numbers returned to get the power |
30 |
;; ratios. |
31 |
|
32 |
;; Because the P and S polarizations are independent, calculate them |
33 |
;; separately. |
34 |
|
35 |
;;; Fresnel's equations from |
36 |
;;; http://chsfpc5.chem.ncsu.edu/CH795Z/lecture/lecture8/fresnel/fresnel.html |
37 |
;;; Altered signs to match |
38 |
;;; http://hyperphysics.phy-astr.gsu.edu/hbase/phyopt/freseq.html |
39 |
|
40 |
;;; Rewrote using matrix method from [Sernelius] |
41 |
;;; http://www.ifm.liu.se/~boser/elma/Lect13.pdf |
42 |
|
43 |
;;; The transmitted voltage (E-Field) |
44 |
(define (E_T n1 n2 cos-i cos-t s-polarization?) |
45 |
(if s-polarization? |
46 |
(/ (* 2 n1 cos-i) |
47 |
(+ (* n1 cos-i) (* n2 cos-t))) |
48 |
(/ (* 2 n1 cos-i) |
49 |
(+ (* n1 cos-t) (* n2 cos-i))))) |
50 |
;;; The reflected voltage (E-Field) |
51 |
(define (E_R n1 n2 cos-i cos-t s-polarization?) |
52 |
(if s-polarization? |
53 |
(/ (- (* n1 cos-i) (* n2 cos-t)) |
54 |
(+ (* n1 cos-i) (* n2 cos-t))) |
55 |
(/ (- (* n2 cos-i) (* n1 cos-t)) |
56 |
(+ (* n2 cos-i) (* n1 cos-t))))) |
57 |
|
58 |
;;; Returns 2x2 matrix: |
59 |
;;; |
60 |
;;; ( 1 r ) |
61 |
;;; 1 ( n-1,n ) |
62 |
;;; -------- * ( ) |
63 |
;;; t ( r 1 ) |
64 |
;;; n-1,n ( n-1,n ) |
65 |
(define (layer-interface n1 n2 th1 s-polarization?) |
66 |
(define cos-i (cos th1)) |
67 |
(define cos-t (cos (Snell-law n1 n2 th1))) |
68 |
(let ((transmit (E_T n1 n2 cos-i cos-t s-polarization?)) |
69 |
(reflect (E_R n1 n2 cos-i cos-t s-polarization?))) |
70 |
(let ((r/t (/ reflect transmit)) |
71 |
(tinv (/ transmit))) |
72 |
(list (list tinv r/t) (list r/t tinv))))) |
73 |
|
74 |
;;; Returns 2x2 matrix coding phase difference between reflected and |
75 |
;;; transmitted paths. |
76 |
;;; |
77 |
;;; ( -i*d_n 0 ) |
78 |
;;; ( e ) |
79 |
;;; ( i*d_n ) |
80 |
;;; ( 0 e ) |
81 |
(define (layer-phase h_j n_j th_j w) |
82 |
(define phase (exp (/ (* +2i pi h_j n_j (cos th_j)) w))) |
83 |
(list (list (/ phase) 0) (list 0 phase))) |
84 |
|
85 |
;;; LAYERS are lists: (index-of-refraction height). The |
86 |
;;; index-of-refraction may be a complex number or a procedure of |
87 |
;;; wavelength (W) returning a complex number. |
88 |
;;; |
89 |
;;; IR_0 is the index-of-refraction of the medium on the top of the |
90 |
;;; stack of LAYERS. The bottom layer has thickness 0. |
91 |
;;; W is probe wavelength. |
92 |
;;; TH_I is the angle of the incident ray. |
93 |
;;; |
94 |
;;; `combine-layers' returns a list of the transmitted and forward and |
95 |
;;; reverse reflected power ratios. The transmitted and each |
96 |
;;; reflected power ratio sum to 1.0 if stack is all dielectric |
97 |
;;; (lossless) |
98 |
;;; |
99 |
;;; Negative W computes the S-polarization, else the P-polarization. |
100 |
(define (combine-layers th_i w layers) |
101 |
(define (numberize x) (if (procedure? x) (x w) x)) |
102 |
(define IR_0 (caar layers)) |
103 |
(define (add-next-0-layer stack) |
104 |
(cons (cons (caadr stack) (cons 0 (cddadr stack))) |
105 |
(cdr stack))) |
106 |
(set! IR_0 (numberize IR_0)) |
107 |
(let ((th_t th_i) |
108 |
(IR_n #f) |
109 |
(M '((1 0) (0 1))) |
110 |
(M2 '((1 0) (0 1))) |
111 |
(s-polarization? (negative? w))) |
112 |
(define (loop stack IR_n-1) |
113 |
(cond |
114 |
((null? stack) |
115 |
(check-finite M2) |
116 |
(matrix->powers (matrix-product |
117 |
M2 (map (lambda (row) (map squmag row)) M)) |
118 |
(magnitude (/ (* (real-part IR_n) (cos th_t)) |
119 |
(real-part IR_0) (cos th_i))))) |
120 |
(else |
121 |
(let ((layer (car stack))) |
122 |
(set! IR_n (numberize (car layer))) |
123 |
(cond |
124 |
((negative? (cadr layer)) ; layer* |
125 |
(set! M2 |
126 |
(matrix-product |
127 |
M2 (map (lambda (row) (map squmag row)) M))) |
128 |
(set! th_t (Snell-law IR_n-1 IR_n th_t)) |
129 |
(set! M2 |
130 |
(matrix-product |
131 |
M2 (map (lambda (row) (map squmag row)) |
132 |
(layer-phase (abs (cadr layer)) IR_n th_t w)))) |
133 |
(set! M '((1 0) (0 1))) |
134 |
(loop (if (negative? (cadadr stack)) ; thickness of next layer |
135 |
;; then add 0-thickness layer with next IR |
136 |
(add-next-0-layer stack) |
137 |
;; otherwise, next layer. |
138 |
(cdr stack)) |
139 |
IR_n)) |
140 |
(else ; layer |
141 |
(set! M (matrix-product |
142 |
M (layer-interface IR_n-1 IR_n th_t s-polarization?))) |
143 |
(set! th_t (Snell-law IR_n-1 IR_n th_t)) |
144 |
(set! M (matrix-product M (layer-phase (cadr layer) IR_n th_t w))) |
145 |
(loop (if (and (positive? (cadr layer)) |
146 |
(negative? (cadadr stack))) |
147 |
(add-next-0-layer stack) |
148 |
(cdr stack)) |
149 |
IR_n))))))) |
150 |
(set! w (abs w)) |
151 |
(loop (if (negative? (cadadr layers)) |
152 |
(add-next-0-layer layers) |
153 |
(cdr layers)) |
154 |
IR_0))) |
155 |
|
156 |
;;; matrix->powers returns a list of the transmitted and forward and |
157 |
;;; reverse reflected power ratios. |
158 |
(define (matrix->powers M IR_ratio) |
159 |
(list (/ IR_ratio (caar M)) |
160 |
(/ (caadr M) (caar M)) |
161 |
(/ (cadar M) (caar M)))) |
162 |
|
163 |
(define (matrix-product mat1 mat2) |
164 |
(map (lambda (arow) |
165 |
(apply map |
166 |
(lambda bcol (apply + (map * bcol arow))) |
167 |
mat2)) |
168 |
mat1)) |
169 |
|
170 |
(define (check-finite M) |
171 |
;;(let ((det (- (* (caar M) (cadadr M)) (* (cadar M) (caadr M))))) ) |
172 |
(if (or (infinite? (caar M)) |
173 |
(infinite? (cadadr M)) |
174 |
(infinite? (cadar M)) |
175 |
(infinite? (caadr M))) |
176 |
(slib:warn 'infinite M)) |
177 |
M) |
178 |
|
179 |
;;(not (and (< .99 (real-part det) 1.01) (< -.01 (imag-part det) .01))) |
180 |
|
181 |
;;; Square of the magnitude |
182 |
(define (squmag x) |
183 |
(define mag (magnitude x)) |
184 |
(* mag mag)) |
185 |
|
186 |
;;; J. C. Maxwell Garnett, |
187 |
;;; "Colours in Metal Glasses and in Metallic Films", |
188 |
;;; Phil. Trans. Roy. Soc. London 203A, 385 (1904).@* |
189 |
;; (define (gran-IR IR q IR0) |
190 |
;; (define e (* IR IR)) |
191 |
;; (define es (* IR0 IR0)) |
192 |
;; (sqrt (* es (/ (+ e (* 2 es) (* -2 (- es e) q)) |
193 |
;; (+ e (* 2 es) (* (- es e) q)))))) |
194 |
(define (gran-IR IR q IR0) |
195 |
(define e (* IR IR)) |
196 |
(define es (* IR0 IR0)) |
197 |
(let ((r (* q (/ (- e es) (+ e es es))))) |
198 |
(sqrt (* es (/ (+ 1 r r) (- 1 r)))))) |
199 |
|
200 |
(define (granular-IR IR q IR0) |
201 |
(if (and (number? IR) (number? IR0)) |
202 |
(gran-IR IR q IR0) |
203 |
(lambda (w) |
204 |
(gran-IR (if (number? IR) IR (IR w)) |
205 |
q |
206 |
(if (number? IR0) IR0 (IR0 w)))))) |
207 |
|
208 |
;; Given the angle of impinging light, |
209 |
;; returns the angle of the transmitted light |
210 |
(define (Snell-law n1 n2 th-i) |
211 |
(asin (* (/ n1 n2) (sin th-i)))) |
212 |
|
213 |
;; Given the angle of impinging light, |
214 |
;; returns the angle of the transmitted light |
215 |
(define (Snells-law n1 n2 th-i) |
216 |
(define sin-th (real-sin th-i)) |
217 |
(define nrat (/ n2 n1)) |
218 |
(real-atan (/ sin-th |
219 |
(real-part (sqrt (- (* nrat nrat) (* sin-th sin-th))))))) |
220 |
|
221 |
(define (find-angles th_0 layers) |
222 |
(do ((layers layers (cdr layers)) |
223 |
(angles (list th_0) |
224 |
(cons (Snells-law (caar layers) (caadr layers) (car angles)) |
225 |
angles))) |
226 |
((null? (cdr layers)) (reverse angles)))) |
227 |
;;(trace-all "fresneleq.scm") (set! *qp-width* 333) |