A function f
f
(x) [ V(
f
) may be represented as
f
f
(x) ¼
X
n[Z
c
n
f
n
(x)(3:124)
In contrast with the sampling theorem and with the Haar wavelet expansion, the
expansion coefficients are not samples of f
f
or inner products between f
f
and the
basis vectors. For the B-splines it turns out that we can derive complementary
functions

f
n
(x) for each
f
n
(x) ¼
b

lementary functions, f
est
is by design consistent with the measurements. We can
further state that f
est
¼ f
f
if the complementary functions are such that

f
[ V(
f
),
in which case there exist discrete coefficients p(k) such that

f
(x) ¼
X
k[Z
p(k)
f
(x Àk)(3:126)
Using the convolution theorem, the Fourier transform of

f
(x)is
^

f
(u) ¼

where a(n) ¼h
f
(x)j
f
(x Àn)i. Without loss of generality, we set n ¼ 0 and sum both
sides of Eqn. (128) against the discrete kernel e
Ài2
p
n
0
u
to obtain
X
n
0
[Z
d
0n
0
e
i2
p
n
0
u
¼ 1
¼
X
k[Z
X

"#
(3:129)
where we use the substitution of variables n
00
¼ n
0
À k.
3.9 B-SPLINES 91
Poisson’s summation formula is helpful in analyzing the sums in Eqn. (3.129).
The summation formula states that for g(x) [ L
1
(R)
X
n[Z
g(n)e
Ài2
p
nu
¼
X
k[Z
^
g(u þk)(3:130)
where
^
g(u) is the Fourier transform of g(x). To prove the summation formula,
note that
h(u) ¼
X
k[Z

ð
kþ1
k
^
g(u)e
2
p
inu
du
¼
ð
1
À1
^
g(u)e
2
p
inu
du
¼ g(n)(3:132)
Since a(x) is the autocorrelation of
f
, its Fourier transform is j
^
f
(u)j
2
. Thus by the
Poisson summation formula
X

1
P
k[Z
j
^
f
(u þk)j
2
(3:134)
92 ANALYSIS
Substitution in Eqn. (3.127) yields
^

f
(u) ¼
^
f
(u)
P
k[Z
j
^
f
(u þk)j
2
(3:135)
We can evaluate Eqn. (3.135) to determine
^

f

2
¼ 1, Eqn. (135)
reduces to
^

f
(u) ¼
^
f
(u) and an orthonormal basis may be obtained.
The Fourier transform of the mth-order B-spline is
^
b
m
(u) ¼ [sinc(u)]
(mþ1)
e
Ài
pj
u
¼
^
b
0
(u)
hi
(mþ1)
e
i
p

jsinc(u þk)j
2
,
meaning that Q
m
(u) Q
o
(u). Thus, 0 , Q
m
(u) , 1 and the B-spline functions of
all orders satisfy the Riesz basis condition.
In contrast with the B-splines themselves, the complementary functions

f
(x)do
not have finite support. It is possible, nevertheless, to estimate

f
(x) over a finite
interval for each B-spline order by numerical methods. Estimation of Q
m
(u)from
Eqn. (3.138) is the first step in numerical analysis. This objective is relatively
easily achieved because Q
m
(u) is periodic with period 1 in u. Evaluation of the
sum over the first several thousand orders for closely spaced values of 0 ! u 1
takes a few seconds on a digital computer.
Given Q
m

f
(x), we can calculate f
f
(x) for target functions. For
example, Fig. 3.17 shows the signals of Figs. 3.8 and 3.9 projected onto the V(
f
) sub-
spaces for B-splines of orders 0–3. Higher-order splines smoothly represent signals
with higher-order local polynomial curvature. Note that higher-order splines are not
more localized than the lower-order functions, however, and thus do not immediately
translate into higher signal resolution. Notice also the errors at the edges of the signal
windows in Fig. 3.17. These arise from the boundary conditions used to truncate the
infinite time signal f (x). In the case of these figures, f (x) was assumed to be periodic
in the window width, such that sampling and interpolation functions extending
beyond the window could be wrapped around the window.
The interpolated signals plotted in Fig. 3.17 are the projections f
f
(x) [ V(
f
)of
f (x) onto the corresponding subspaces V(
f
). The consistency requirement designed
into the interpolation strategy means that these functions, despite their obvious discre-
pancies relative to the actual signals, would yield the same sample projections.
Corrections that map the interpolated signals back onto the actual signal lie in
V
?
(
f


f
. As the order of the
B-spline tends to infinity,

f
(x) converges on sinc(x) [235]. If we remove the require-
ment that

f
(x) [ V(
f
), it is possible to generate a biorthogonal dual basis for
b
m
(x)
with compact support [49]. The compactly supported biorthogonal wavelets in this
case introduce a complementary subspace

V spanned by

f
(x).
The goal of the current section has been to consider how one might use a set of
discrete B-spline inner products to estimate a continuous signal. This problem is
central to imaging and optical signal analysis. We have already encountered it in
the coded aperture and tomographic systems considered in Chapter 2, and we will
encounter it again in the remaining chapters of the text. We leave this problem for
now, however, to consider the use of sampling functions and multiscale represen-
tations in signal and system analysis. One may increase the resolution and fidelity

Wavelet theory is a broad and powerful branch of mathematics, and the student is
well advised to consult standard courses and texts for deeper understanding [53,164].
Wavelets often describe images and other natural signals well. The intuitive match
between wavelets and images arises from the assumption that “features” in natural
signals tend to cluster, meaning that higher resolution is desirable in the vicinity of
a feature than elsewhere in the signal. Multiscale clustering enables wavelet represen-
tations to estimate signals with fewer samples than might be used with uniform
regular sampling. Under the Whittaker–Shannon sampling strategy, functional
samples are distributed uniformly in space even in regions with no significant
image features. Wavelets enable samples to be dynamically assigned to regions
with interesting features. This dynamic resource allocation is the basis of natural
signal compression.
B-splines may be used to generate semiorthogonal bases as in Section 3.9, biortho-
gonal spaces and orthogonal wavelet bases. As before, we imagine a hierarchy
of spaces
{0} , ÁÁÁ, V
2
, V
1
, V
o
, V
À1
, V
À2
, ÁÁÁ, L
2
(R)(3:140)
Semiorthogonal bases are spanned by sets of functions that are not themselves
orthogonal but are orthogonal to a complementary set of functions. Biorthogonal

^
p(u)
^
b
m
(u)(3:142)
Our goal is to select
f
(x) to be an orthonormal scaling function such that
f
(x Àn),
f
(x À m)
hi
¼
ð
1
À1
f
Ã
(x Àn)
f
(x Àm)dx
¼
d
nm
(3:143)
We may apply the Poisson summation formula as in Section 3.9 to derive a simple
identity from Eqn. (3.143). Again letting a(x) ¼
f

a(u) ¼j
^
f
(u)j
2
, which yields the identity for orthonormal scaling functions
X
k
j
^
f
(u þ k)j
2
¼ 1(3:145)
Referring to Eqn. (142), we see that
^
f
(u) satisfies Eqn. (145) if we select
^
p(u) ¼
1
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
P
k
j
^
b
m
(u þk)j
2

where
S
n
(u) ¼
X
k[Z
1
(u þ k)
n
(3:148)
We know that the m ¼ 0 spline produces the Haar scaling function
^
f
0
(u) ¼ e
Ài
p
u
sin (
p
u)
p
u
(3:149)
Comparing Eqns. (3.147) and (3.149), we see that
S
2
(u) ¼
p
2

(u) ¼
p
6
(33 þ26 cos (2
p
u) þcos (4
p
u))
180 sin
6
(
p
u)
(3:152)
and
S
8
(u) ¼
p
8
(1208 þ1191 cos (2
p
u) þ120 cos (4
p
u) þ cos(6
p
u))
10,080 sin
8
(

ffiffiffiffiffiffiffiffiffi
2
jÀ1
p
h[n
0
À n]
f
x
2
jÀ1
À n
0

(3:154)
Equation (3.154) reduces without loss of generality to
1
ffiffiffi
2
p
f
x
2
¼
X
n
h[n]
f
x À nðÞ (3:155)
The Fourier transform of Eqn. (155) yields

(2u)
^
f
(u)
¼ e
Ài
pj
u
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
S
2 mþ2
(u)
2
2 mþ1
S
2 mþ2
(2u)
s
(3:158)
As with the Haar scaling function, we are interested in obtaining orthogonal wave-
lets spanning the spaces W
j
such that V
jÀ1
¼ V
j
È W
j
. Such wavelets are immediately
obtained using the conjugate mirror filter

orthogonality and scaling rules as the Haar wavelets discussed earlier; specifically
f
j,n
(x) ¼
1
ffiffiffiffiffi
2
j
p
f
x
2
j
À n

(3:160)
c
j,n
(x) ¼
1
ffiffiffiffiffi
2
j
p
c
x
2
j
À n


0
d
nn
0
(3:163)
c
j,n
j
f
j
0
,n
0

¼ 0(3:164)
As with the Haar wavelets, the Battle–Lemarie functions span L
2
(R
2
) in the hier-
archy of spaces described by Eqn. (3.140). The Battle–Lemarie wavelets are pre-
sented here to provide an accessible introduction to wavelet theory. Many other
wavelet families have been developed [164]; the selection of which family to use
for a particular class of signals is application-specific. Some wavelets are attractive
because they have compact support, which the Shannon wavelet famously does
not. Other wavelets, such as the Haar and B-splines, arise naturally from physical
3.10 WAVELETS 99
or system design considerations. In still other cases, a particular basis may prove more
amenable to compact support of a particular signal class.
PROBLEMS

e
x
2
d
n
dx
n
e
Àx
2
(3:166)
Defining
f
n
(x) ¼ e
À
p
x
2
H
n
(
ffiffiffiffiffiffi
2
p
p
x)(3:167)
show that for n . 0
ffiffiffiffiffiffi
2

Can you describe the structure of the aliased signal?
(b) Plot the discrete Fourier transform of sin (2
p
ux) on [0, 1] using 1024 uni-
formly spaced samples for u ¼ 16, 32, 64. Label the plot in frequency
units. What is the width of the Fourier features that you observe? What
causes this width?
(c) Plot the discrete Fourier transform of
b
0
(x) sin (2
p
ux)on[À1:5, 2:5]
using 4096 uniformly spaced samples for u ¼ 16, 32, 64. Label the plot
in frequency units and explain the plot.
100 ANALYSIS
3.6 Fourier Analysis of a Hermite–Gaussian:
(a) Plot the Hermite–Gaussian
f
5
(x) over the range of the function.
(b) Plot
~
f
5
t
(x) over a representative range of
t
.
3.7 Transformations:

(b) Replicate Figs. 3.9, 3.12, and 3.17 for your function.
3.11 2D Wavelet Analysis. Replicate Figs. 3.13 and 3.14 for an image of your
choosing.
3.12 Spline Interpolation:
(a) Show that a one-dimensional pinhole imaging system produces
measurements
g
n
¼ f
n
¼
ð
f (x)
b
1
x À nD
D

dx (3:170)
(b) Plot the values of f
n
for
f (x) ¼ cos 2
p
x
5D

e
À(x=30D)
2

coherence properties of the source are central to our discussion. Coherence theory,
which relates the electromagnetic nature of the field to statistical properties of
quantum (e.g., photonic) processes, is the focus of Chapter 6.
†
How the field is detected. The field may be detected by optically induced
chemical, physical, thermal, and electronic effects. Optoelectronic detection
interfaces for imaging and spectroscopy are the focus of Chapter 5.
Optical Imaging and Spectroscopy. By David J. Brady
Copyright # 2009 John Wiley & Sons, Inc.
103
†
How the field propagates and how propagating fields are modulated by
materials. Field propagation is described by the Maxwell equations for electro-
magnetic waves, and field–matter interactions are described by materials
equations. The electromagnetic description of optical waves and optical inter-
actions is the focus of this chapter.
In view of the peculiarly quantum mechanical nature of optical field generation and
detection, it is important to understand that the conventional electromagnetic field of
the Maxwell equations is not a sufficient description of optical fields. The description
derived in this chapter provides a basis for optical analysis, but complete understand-
ing of optical field propagation and field properties must incorporate the detection
and coherence processes discussed in Chapters 5 and 6. In short, the student must
understand the next three chapters as a group to have a vision for the peculiar and
beautiful nature of optical fields.
4.2 WAVE MODEL FOR OPTICAL FIELDS
The Maxwell equations for electromagnetic propagation are
rÂE ¼À
@
@t
B (4:1)

linear dielectric relationship
P ¼ 1
0
x
e
E (4:7)
such that
D ¼ 1E (4:8)
where
1 ¼ (1 þ
x
e
)1
0
(4:9)
Since charge dynamics at optical frequencies are described by quantum mechanical
processes that cannot be accurately analyzed by continuous models, the space charge
density
r
and the current density J are generally neglected in optical analysis. A
nonzero current density is sometimes applied to formally account for optical absorp-
tion. We also note that E and D need not be collinear, meaning that 1 is in general
tensor-valued. Materials in which 1 is a scalar are called isotropic. Optical glasses
are isotropic, but optical crystals are often anisotropic. While most of the optical
systems discussed in this text utilize isotropic materials, we consider the application
of anisotropic materials to tunable filters in Section 9.7.
Using the material relations to eliminate B and D from the Faraday and Ampere
relationships yields
rÂE ¼À
m

A (4:14)
From Gauss’ law we know that
rÁ 1EðÞ¼E Ár1 þ 1rÁE ¼ 0(4:15)
where we have assumed for the moment that 1 is a scalar. We can reexpress
Eqn. (4.15) as
rÁE ¼ÀE Árlog 1ðÞ (4:16)
4.2 WAVE MODEL FOR OPTICAL FIELDS 105
Substituting Eqn. (4.16) in the wave equation yields
r
2
E ÀrE Árlog 1
ðÞðÞ
¼
m
0
1
@
2
@t
2
E (4:17)
A medium in which r1 ¼ 0 is termed homogeneous. Most optical dielectrics (i.e.,
transparent glasses and crystals and liquids) are homogeneous. Optically interesting
inhomogeneous media include photonic crystals, optical fiber, graded-index lenses,
and volume holograms.
In isotropic homogeneous media, the wave equations reduce to
r
2
E À
m

u  u  E
0
¼À
m
0
1n
2
E
0
(4:21)
If 1 is a scalar, this equation has solutions E
0
only if juj
2
¼
m
0
1n
2
. Allowing for the
possibility that 1 is a tensor, solutions correspond to values of u such that
u Âu Âþ
m
0
1n
2





expressed as superpositions of plane waves. In the remainder of this chapter we
restrict our attention to isotropic materials, in which case 1 is a scalar and juj
2
¼
m
0
1n
2
. A general solution to Eqn. (4.18) in this case is
E(r, t) ¼
ððð
F(u, v, n)p(u, v, n)e
Ài2
p
ntÀuxÀvyÀ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
m
1n
2
Àu
2
Àv
2
z
p

du dv dn (4:23)
where (u, v, w ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
m

i2
p
uxþvyþ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
m
1n
2
Àu
2
Àv
2
z
p

du dv (4:25)
In the remainder of this chapter we also assume that the field propagates paraxially.
The paraxial approximation consists of the assumption that values of u and v for
which jF(u, v)j is nonzero lie on a compact window on the wave normal sphere
Figure 4.1 The wave normal surface in free space.
4.3 WAVE PROPAGATION 107
centered on the w axis, as illustrated in Fig. 4.1. This window is centered on the z axis,
such that w ) u, v over the full spatial bandwidth. This means that the polarization
vector p(u, v) is nearly parallel to the (x, y) plane over the entire spatial bandwidth.
In an isotropic material p(u, v) may be represented on any basis orthogonal to u.
We select as an example a basis in which one of the polarization vectors is also
orthogonal to the y axis. The resulting orthonormal basis for p(u, v)is
p
x
¼À
k

þ
l
vi
y
þ
l
wi
z

 p
x
¼
(1 À
l
2
v
2
)i
y
þ
l
2
uvi
x
þ
l
2
i
z
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1=
l
2
Àu
2
Àv
2
z
p

du dv
f
y
(r) ¼
ðð
F
y
(u, v)p
y
e
i2
p
uxþvyþ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1=
l
2
Àu
2

x
(x, y, z)e
Ài2
p
(uxþvy)
dx dy
F
y
(u, v)p
y
e
i2
p
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1=
l
2
Àu
2
Àv
2
z
p
¼
ðð
f
y
(x, y, z)e
Ài2
p

y
¼
ðð
f
y
(x, y, z ¼ 0)e
Ài2
p
(uxþvy)
dx dy
(4:29)
Once we have determined F
x
(u, v) and F
y
(u, v), the field at all (x, y, z)maybe
calculated from Eqn. (4.27). Specification of the field on a surface, such as the
108 WAVE IMAGING
plane z ¼ 0, is called a boundary condition and the evolution of the field distribution
from one boundary to another is called diffraction. Equations (4.27) and (4.28) enable
us to computationally model diffraction in homogeneous media.
4.4 DIFFRACTION
Diffraction is the process of wave propagation from one boundary to another. A
canonical example of optical diffraction, propagation of a monochromatic field
from the plane (x, y, z ¼ 0) to the plane (x
0
, y
0
, z ¼ d), is illustrated in Fig. 4.2.
Given the electric field distribution on the input plane, we seek to estimate the

u. The F
z
component
does not produce a propagating field.
Let G
x
(u, v) and G
y
(u, v) be the Fourier distributions of the field in the x
0
, y
0
plane
at z ¼ d. From Eqn. (4.28) we see that
G
x
(u, v) ¼ e
i2
p
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1=
l
2
Àu
2
Àv
2
d
p
F

Àu
2
Àv
2
d
p
is the transfer function for diffraction from the
z ¼ 0 plane to the z ¼ d plane.
Nominally, the impulse response for diffraction is the inverse Fourier transform of
the transfer function. We continue along this line with care, however, by briefly
accounting for the vector nature of the field. Using the transfer function and
Eqn. (4.27), we obtain
g(x
0
, y
0
) ¼
ðð
F(u, v) Áp
x
p
x
þ F(u, v) Áp
y
p
y

 e
i2
p

In this case Eqn. (4.32) simplifies considerably to yield
g(x
0
, y
0
) ¼
ðð
F(u, v)e
i2
p
ux
0
þvy
0
þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1=
l
2
Àu
2
Àv
2
z
p

du dv
¼
ðð
f(x, y)h(x

þ x
2
þ y
2
)
þ
l
2
p
(d
2
þ x
2
þ y
2
)
3=2

e
i(2
p
=
l
)
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
d
2
þx
2
þy

methods described in Section 3.7. In analytic work, however, one generally
110 WAVE IMAGING
chooses to work with a simplified approximate impulse response. As an example,
l
( d in essentially all optical systems, meaning that the 1/
l
d term in Eqn. (4.36)
dominates the 1/d
2
term. In imaging system analysis, the impulse response is
often simplified by the Fresnel (near-field) approximation or the more restrictive
Fraunhofer (far-field) approximation. Both approximations are paraxial, meaning
that we restrict our attention to field distributions over the space close to the axis
of optical propagation (the z axis in Fig. 4.2).
The Fresnel approximation is just the paraxial approximation that d )jx Àx
0
j,
jy À y
0
j for all x, y and x
0
, y
0
of interest. In this case
h(x, y) %
1
i
l
d
e

/
D
. 1, one might expect features
of size
D
to blur on propagation at distances greater than d ¼
D
2
/
l
. This suggests
that wavelength scale features will blur quite rapidly on diffraction. Features with
an initial scale of 10 wavelengths blur in 100 wavelengths, while features on a
scale of 100 wavelengths blur in 10,000 wavelengths. This effect is illustrated
in Fig. 4.3, which shows diffraction of Gaussian spots of various sizes. Notice,
however, that high-frequency features reappear at 10 mm as a result of interference
between the diffracting spots. Such interference appears in the diffraction of coherent
laser fields, but is not observed in the diffraction of incoherent fields.
Figure 4.3 was generated using numerical analysis in Matlab. The figure used a
2 Â 2-mm spatial window sampled with 1024 Â 1024 pixels. The Fresnel transfer
function multiplied the DFT of the input field and an inverse DFT was used to gen-
erate the diffracted field. Of course, numerical analysis is not necessary for analysis of
diffraction of these particular sources because, as discussed in Section 3.5, Hermite–
Gaussian distributions are eigenfunctions of the Fresnel transform. According to
Eqns. (3.76) and (4.38), if the input field f (x, y) ¼
f
n
(x=w
0
)

þ
l
2
d
2
q
e
i
p
½(x
2
þy
2
)
l
d=(w
4
0
þ
l
2
d
2
)
Â
f
n
xw
0
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

B
@
1
C
A
(4:39)
With specific reference to the Gaussian input spots of Fig. 4.3, the diffracted field for
the input distribution f(x, y) ¼ exp (À
p
(x
2
þ y
2
)=w
2
0
)is
g(x, y) ¼
e
Ài(
p
=4)
w
2
0
þ i
l
d
e
i(2

to note that blur is not a fundamental process of diffraction for coherent fields. In fact,
a diffracting coherent field maintains its spatial frequency bandwidth on propagation.
Figure 4.3 Absolute values of the diffracted field for the Gaussian spots e
À
p
½(x
2
þy
2
)=20
2

,
e
À
p
f½(xÀ100)
2
þy
2
=50
2
g
, e
À
p
f½(xþ200)
2
þy
2

e
i(2
p
d=
l
)
i
l
d
ðð
e
i(
p
=
l
d)
[(x À x
0
)
2
þ ( y À y
0
)
2
] f (x, y)dx dy
¼
e
i(2
p
d=

(x
2
þ y
2
)
hi
f (x, y)dx dy (4:41)
Figure 4.4 Absolute value of the diffracted field for f (x, y) ¼ e
À
p
[(x
2
þy
2
)=250
2

½1 þ
cos (0:05
p
x)]. All units are in micrometers;
l
¼ 1 mm.
4.4 DIFFRACTION 113
Assuming that x
2
(
l
d and y
2

^
fu¼
x
0
l
d
, v ¼
y
0
l
d

(4:42)
meaning that the diffracted field is proportional to the Fourier transform of the input
field. Note that the Fraunhofer assumption is quite restrictive, however. For example,
a 100 wavelength scale input must diffract for well over 10,000 wavelengths to reach
the Fraunhofer regime and a 1000 wavelength feature, for well over 1,000,000 wave-
lengths. Fraunhofer diffraction is nevertheless often useful in determining the rough
size and spatial frequency structure of objects. The Fraunhofer assumption is com-
monly applied at opposite ends of the electromagnetic imaging frequency scale,
such as in X-ray crystallography and radio astronomy.
Figure 4.5 Absolute value of the diffracted field for f (x, y) ¼ e
À
p
½(x
2
þy
2
)=250
2

(r) ¼ E
i
exp(2
p
iu
i
.
r).
The incident wave is refracted at the interface into the second medium of index of
refraction n
2
, and a reflected wave is returned into the first medium. The refracted
and reflected fields are E
t
(r) ¼ E
t
exp(2
p
iu
t
.
r) and E
r
(r) ¼ E
r
exp(2
p
iu
r
.

r
)] Âi
s
¼
E
t
(
r
) Â i
s
, for all
r
on the interface. i
s
is the surface normal for the interface. To
satisfy this condition, one must require that
u
i
À u
i
Á i
s
i
s
¼ u
t
À u
t
Á i
s

2
1
l
2
þ u
i
Á i
s
ðÞ
2
s
(4:44)
If we use angles relative to the surface normal to decompose u
i
and u
t
into trans-
verse and longitudinal components, Eqn. (4.44) immediately reduces to Snell’s law
[Eqn. (2.5)]. In the the paraxial case i
s
and u
i
are nearly collinear and Eqns. (4.44)
can be approximated by
u
t
% u
i
þ i
s

is the surface normal at the first interface of a prism and i
2
the surface
normal at the second interface, recursive application Eqn. (4.45) produces an estimate
of the output wavevector u
3
u
3
% u
2
À i
2
Dn
"
n
l
2
u
2
Á i
2
% u
1
þ
Dn
"
n
l
(i
1