/[freesnell]/freesnell/fresneleq.scm
ViewVC logotype

Contents of /freesnell/fresneleq.scm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.20 - (show annotations) (download)
Fri Mar 23 04:18:29 2012 UTC (12 years, 1 month ago) by jaffer
Branch: MAIN
CVS Tags: FreeSnell-1c1, FreeSnell-1c2, FreeSnell-1c3, FreeSnell-1c4, FreeSnell-1b3, HEAD
Changes since 1.19: +2 -2 lines
From Vitali Mueller
(combine-layers): Corrected transmission for non-normal incidence with
unequal substrate IRs.

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)

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26