Ersoy O.K. Diffraction, Fourier optics, and imaging (Wiley, 2006)(ISBN 0471238163)(427s) PEo

88	INVERSE DIFFRACTION

Depending on the values of M and Q, there are two possible cases that are discussed below.

Case I: Q M

In this case, Eq. (6.4-5) reduces to

Uðx; zrÞ ¼ 1	2 M	ej2pfx0ðx x0Þdfx03Tðx; zrÞdx
ð	ð
1	4 M			5
¼ Tðx; zrÞ sin c			2Mx
			l
¼ Tðx; zrÞ
because T(x,zr) is bandlimited.
Case II: Q < M
In this case, Eq. (6.4-12) reduces to
Uðx; zrÞ ¼ Tðx; zrÞ sin c					Qx
					2
					l

ð6:4-5Þ

ð6:4-6Þ

This means lowpass filtering of T(x,zr) because the Fourier transform of the sin c function is the rectangular function. In other words, frequencies above Q=l are filtered out. The reconstructed image may have smoothed edges and loss of detail as a consequence.

As Q is restricted to the values Q 1=l in practice, the resolution achievable in the image is l. This is in agreement with other classical results that indicate that linear imaging systems cannot resolve distances less than a wavelength.

EXAMPLE 6.1 Determine

U(x,zr)

when

Uðx0; z0Þ is measured in

an aperture

jx0j R.

Solution: In this case, Eq. (6.4-3) becomes

ð

x; zr

Þ ¼

ð

8e jkz0r pk2 4p2fx2

ð

dx0e j2pfxx0

U

Q

R

Q

<

R

dxTðx; zrÞe j2pfx0x3319

02

dfx0ej2pfx0x0 2

1

dfx

:

M

ð

=

j2pfxx

@4

ð

4

55A

;

M 1

ANALYSIS

89

Interchanging orders of integration, this can be written as
Uðx; zrÞ ¼	ð		ð		dfx0ðe j2pfx0xejkz0r pk2 4p2fx02
	1	2 M
	1	4		M
				Q		ej2pfxxe jz0r pk2 4p2fx2 2R sin c 2Rfx		3T x; zr dx
					dfx
				Q			Þ&Þ	7	ð	Þ
				ð	½	ð		5

It is observed that the sin c function above mixes the plane wave components in the original object wave. This means that each plane wave component in the original object wave is truncated by the aperture (-R,R) and therefore ‘‘spreads’’ by diffraction.

EXAMPLE 6.2 When M ¼ Q ¼ 1=l, show that T ðx; zrÞ can be exactly recovered from Uðx0; z0Þ by forward propagation of U ðx0;z0Þ instead of inverse propagation. Solution: Suppose that forward propagation from the plane z ¼ z0 to z ¼ 2z0 zr is carried out. As

2z0 zr z0 ¼ z0r

the forward propagation equation is given by

				Q		1
			Q			1
Uðx; 2z0 zrÞ ¼				ð dfx		ð		U ðx0; z0Þe j2pfxx0 ejz0r pk2 4p2fx2 ej2pfxxdx0
Rearranging orders of integration results in
U x; 2z0	zr	2 Q	1 U x0; z0					ej2pfxx0 e jz0r pk2 4p2fx2 e j2pfxxdx03 dfx T			x; 0
ð	Þ ¼	6 Q 1			ð		Þ		7	¼ ð Þ
		ð	ð						5
		4

In other words, the original image can be reconstructed by back space propagation of the recorded wave or forward propagation of the complex conjugate of the recorded wave.

7

Wide-Angle Near and Far Field Approximations for Scalar Diffraction

7.1INTRODUCTION

Approximations for computing forward and inverse diffraction integrals are of vital significance in many areas involving wave propagation. As discussed in Chapters 4 and 5, approximations such as the Fresnel approximation, the Fraunhofer approximation, and the more rigorous angular spectrum method (ASM) all involve the Fourier transform, its discrete counterpart, the discrete Fourier transform, and its fast computational routine, the fast Fourier transform (FFT).

The Fresnel approximation is valid at reasonable distances from the input plane whereas the Fraunhofer approximation is valid in the far field. The ASM is a rigorous solution of the Helmholtz equation; its numerical implementation is usually done with the FFT, and possibly other digital signal processing algorithms, with their related approximations [Mellin and Nordin, 2001; Shen and Wang, 2006].

When the sizes of the diffracting apertures are less than the wavelength, scalar diffraction theory yields nonnegligible errors, and other numerical methods such as the finite difference time domain (FDTD) method and the finite element method (FEM) may become necessary to use [Kunz, 1993; Taflove and Hagness, 2005]. However, these methods are not practical with large scale simulations as compared with methods utilizing the FFT. With diffracting aperture sizes of the order of the wavelength used and in the near field, the ASM has been found to give satisfactory results [Mellin and Nordin, 2001]. With the ASM, one disadvantage is that the input and output plane sizes are the same. The output plane size is usually desired to be considerably larger than the input plane size in most applications, and this can be done with the ASM only by repeatedly using the ASM in short distances and additionally using filtering schemes to make the output size progressively larger.

The Fresnel, Fraunhofer approximations, and the ASM can be considered to be valid in practice at small angles of diffraction. It is desirable to have approximate methods that are valid at wide angles of diffraction in the near field as well as the far field and also based on the Fourier integral to be implemented with the FFT so that

Diffraction, Fourier Optics and Imaging, by Okan K. Ersoy

Copyright # 2007 John Wiley & Sons, Inc.

90

A REVIEW OF FRESNEL AND FRAUNHOFER APPROXIMATIONS

91

large scale computations can be carried out in a reasonable amount of time and storage. Even if the spatial frequencies are sampled such that the FFT cannot be used, the Fourier integral representation is still desirable as it is separable in the two variables of integration, reducing 4-D tensor operations to 2-D matrix operations.

In this chapter, a new set of approximations having such features is discussed. The approximations are based on the Taylor expansion around the radial distance of a point on the output plane from the origin. With these approximations, what makes possible to utilize the FFT is semi-irregular sampling at the output plane. Semiirregular sampling means sampling is first done in a regular array, and the array points are subsequently perturbed to satisfy the FFT conditions.

Some interesting results will be related to the similarity of the proposed approximation to the Fresnel approximation even though the new approximation is valid in the near field and at wide angles whereas the Fresnel approximation is not. This is believed to be one reason why the Fresnel approximation has often been used in the near field in many applications. The results obtained in such studies are actually valid at the output sampling points perturbed in position, as discussed in this paper.

This chapter consists of eight sections. Section 7.2 is a review of Fresnel and Fraunhofer approximations previously discussed in Chapter 5. Section 7.3 introduces the new methods with the radial set of approximations. Section 7.4 provides further higher order improvements and error analysis. Section 7.5 shows how the method can be used in inverse diffraction and iterative optimization applications. Section 7.6 provides numerical simulation examples in 2-D and 3-D geometries. Section 7.7 is about how to increase accuracy by centering input and output plane apertures around some center coordinates and possibly also using smaller subareas in large scale simulations. Section 7.8 covers conclusions.

7.2 A REVIEW OF FRESNEL AND FRAUNHOFER APPROXIMATIONS

The first Rayleigh–Sommerfeld diffraction integral can be written as

Uðx; y; zÞ ¼ j1l		ðð	Uðx; y; 0Þ rz2 e jkrdxdy			ð7:2-1Þ
		1
		1

where Uðx; y; 0Þ is the input field, Uðx0; y0; zÞ is the output field, k is the wave number equal to 2p=l, l is the wavelength, and r is the radius vector length from the point ðx; y; 0Þ to the point ðx0; y0; zÞ. In general, z is often large enough to approximate z=r2 by 1=z. Then, Eq. (7.2-1) becomes

			Uðx; y; zÞ ¼ jlz	ðð	Uðx; y; 0Þe jkrdxdy	ð7:2-2Þ
			1	1
				1
The constant term	1		will be neglected in the rest of the chapter.
	jlz

92	WIDE-ANGLE NEAR AND FAR FIELD APPROXIMATIONS

In the Fresnel approximation commonly used in diffraction problems, the radius vector length is approximated by the Taylor’s series expansion with two terms [Mezouari and Harvery, 2003; Southwell, 1981; Steane and Rutt, 1989]:

q					þ	g		ð		Þ
r ¼ z2 þ ðx0 xÞ2 þ ðy0 yÞ2				z 1		2			7:2-3
where
g	¼	ðx0 xÞ2 þ ðy0 yÞ2						ð	7:2-4	Þ
		z2

The magnitude of the maximum phase error Emax in radians is estimated by the next term in the Taylor series as

Emax	1	ð7:2-5Þ
	8 zg2k

The Fresnel approximation is not sufficiently good in problems where Emax is larger than, say, 1 radian. However, this condition has often been relaxed in various applications because good results have been generally obtained. This was partially explained by using the method of stationary phase [Goodman]. In this chapter, another explanation for these results is provided.

The advantage in the Fresnel approximation is that the variables are separable so that the integral in Eq. (7.2-2) can be written as a Fourier integral. For example, Eq. (7.2-4) is rewritten as

g	¼	v w				2ðx0 x þ y0yÞ
		z2	þ	z2		z2

where

v ¼ x20 þ y20

w ¼ x2 þ y2

ð7:2-6Þ

ð7:2-7Þ

ð7:2-8Þ

Then, with the help of Eq. (7.2-3), Eq. (7.2-2) can be approximated as a Fourier integral by

ðð

Uðx0; y0; zÞ ¼ e jk2vz u0ðx; yÞe jkz ðx0xþy0yÞdxdy ð7:2-9Þ

where

u0ðx; yÞ ¼ uðx; yÞe jk	w	ð7:2-10Þ
	2z

THE RADIAL SET OF APPROXIMATIONS

93

Equation (7.2-9) becomes a Fourier integral by writing it as

ðð

Uðx0; y0; zÞ ¼ e jk2vz u0ðx; yÞe j2pðfxxþfyyÞdxdy

where the spatial frequencies fx and fy are defined by

fx ¼ xs lz

fy ¼ ys lz

where xs equals x0, and ys equals y0 at this point.

In the Fraunhofer approximation, the condition z w allows u0ðx; yÞ

that Eq. (7.2-9) becomes
Uðx0; y0; zÞ ¼ e jk2vz	ðð uðx; yÞe jkz ðx0xþy0yÞdxdy

ð7:2-11Þ

ð7:2-12Þ

ð7:2-13Þ

uðx; yÞ so

ð7:2-14Þ

Equation (7.2-9) can be evaluated with the fast Fourier transform [Brigham, 1974] by letting

xs ¼ mx 0x; ys ¼ my 0y;

x¼ nx x;

y¼ ny y;

lz

0x x ¼ Mx

0y y ¼ lz My

mx ¼ integer	ð7:2-15Þ
my ¼ integer	ð7:2-16Þ
nx ¼ integer	ð7:2-17Þ
ny ¼ integer	ð7:2-18Þ
	ð7:2-19Þ
	ð7:2-20Þ

where Mx and My are the number of samples in the x- and y-directions, respectively. In any other approximation, it is necessary to keep Eqs. (7.2-19) and (7.2-20) valid if FFT is to be used.

7.3THE RADIAL SET OF APPROXIMATIONS

In many applications such as those involving digital holography and diffractive optical elements, the coordinates ðx; y; 0Þ are usually much smaller than the object coordinates ðx0; y0; zÞ. Then, Eq. (7.2-3) will be written as

	¼		ð		Þ
r		r0p1 þ h		7:3-1

94 WIDE-ANGLE NEAR AND FAR FIELD APPROXIMATIONS

	where	¼ q				ð		Þ
	r0			z2 þ x02 þ y02			7:3-2
	h		w		2ðx0x þ y0yÞ		7:3-3
		¼ r02			r02	ð		Þ

In order to be able to use the Fourier transform, Eq. (7.3-3) is to be approximated by

	h		w		2ðx0x þ y0yÞ	7:3-4
		¼ z2			r02		Þ
						ð

The simplest radial set of approximations is obtained by approximating the radial distance r as

r		r		w2		xx0 þ yy0	7:3-5
			0 þ	2z				Þ
						r0	ð

Equations (7.3-4) and (7.3-5) are also valid when x; y; x0, and y0 are sufficiently smaller than z (sufficiency conditions will be more clear later). Equation (7.3-5) is further improved in the next section. In order to obtain a Fourier transform relationship with the diffraction integral given by Eq. (7.2-9 ), the following is defined:

x0 ¼	xs	r0	ð7:3-6Þ
	z
y0 ¼	ys	r0	ð7:3-7Þ
	z

where ðxs; ysÞ is the regular output sampling point. The simultaneous solution of Eqs. (7.3-6) and (7.3-7) is given by

x0 ¼ Cðxs; ysÞxs

ð7:3-8Þ

y0 ¼ Cðxs; ysÞys

ð7:3-9Þ

where

C

x

; y

sÞ ¼

1

xs2

þ ys2

21

ð

7:3-10

Þ

z2

ð

s

2

z

w as in the

Fraunhofer

approximation, the

term involving

When

w in Eq. (7.3-5) is skipped. The resulting approximations will be referred to as the near field approximation (NFA) and the far field approximation (FFA) below. The general case covering both NFA and FFA will be referred to as NFFA.

HIGHER ORDER IMPROVEMENTS AND ANALYSIS		95
With the new approximations, Eq. (7.2-11) is replaced by
Uðx0; y0; zÞ ¼ e jkr0	ðð u0ðx; yÞe j2pðfxxþfyyÞdxdy	ð7:3-11Þ

where the spatial frequencies are still given by Eqs. (7.2-12) and (7.2-13) due to Eqs. (7.3-6) and (7.3-7). With FFA, u0ðx; yÞ equals uðx; yÞ.

We note that Eq. (7.3-11) looks like Eq. (7.2-11) for the Fresnel approximation, but the phase factor outside the integral and the output sampling points are different. Equation (31) is straightforward to compute with the FFT, and the results are valid at the sampling points defined by Eqs. (7.3-8) and (7.3-9). The spatial frequencies can also be written as

fx ¼	x0	ð7:3-12Þ
	lr0
fy ¼	y0	ð7:3-13Þ
	lr0

It is observed that fx and fy are still given by Eqs. (7.2-12) and (7.2-13) with respect to the regular sampling points ðxs; ysÞ or by Eqs. (7.3-12) and (7.3-13) with respect to the actual sampling points ðx0; y0Þ.

7.4HIGHER ORDER IMPROVEMENTS AND ANALYSIS

We write r in Eq. (7.3-1) in a Taylor series as

r

¼

r

1

h

h2

h3

1

n 1

1

3 5

ð2n 1Þ

h2m 1

þ

2

8 þ

16 þð

Þ

0

2 4 6 2n

ð7:4-1Þ

where it is desirable to attain the terms shown in order to obtain high accuracy. In order to be able to use the Fourier transform, we only keep the terms in the Taylor series that do not have factors ðx0xÞn other than x0x.

We also define

a	n ¼ ð	1	Þ	n 1	1 3 5 ð2n 1Þ								1
				2			4		6		2n z2n

Then, r is approximated by

r r0 xx0 þ yy0 þ pðwÞ r0

ð7:4-2Þ

ð7:4-3Þ

96 WIDE-ANGLE NEAR AND FAR FIELD APPROXIMATIONS

where

X	an
1
pðwÞ ¼ pðx; yÞ ¼	z2n 1	wn	ð7:4-4Þ
n¼1

Equations (7.3-8) and (7.3-9) remain the same because the term involving xx0 and yy0 is not modified. When z w as in the Fraunhofer approximation, the term pðwÞ in Eq. (7.4-3) is skipped. With the new approximations, Eq. (7.3-11) is replaced by

Uðx0; y0; zÞ ¼ e jkr0	ðð u00ðx; yÞe j2pðfxxþfyyÞdxdy	ð7:4-5Þ
where
u00ðx; yÞ ¼ uðx; yÞe jkpðwÞ		ð7:4-6Þ

For a given z, the maximum error in the Fresnel approximation is proportional to ðx20 þ y20Þ2 whereas, in the NFFA, it is approximately given by

E	¼	vðx2 þ y2Þ		ðx0x þ y0yÞ2	7:4-7	Þ
		4z3		2z3	ð

This error corresponds to the neglected terms in the second component u22 of the Taylor series expansion.

We note that changes are made only in the object coordinates with Eqs. (7.3-8) and (7.3-9). Otherwise, the integral in Eq. (7.4-5) is the same as the integral in Eq. (7.2-9) or Eq. (7.2-14). The output points ðx0; y0Þ are determined by Eqs. (7.3-8) and (7.3-9) in terms of ðxs; ysÞ given by Eqs. (7.2-15) and (7.2-16).

An interesting observation is that the Fresnel approximation examples cited in the literature are approximately correct in the near field at the points ðx0; y0Þ given by Eqs. (7.3-8) and (7.3-9), and not at the points ðxs; ysÞ. As such, these results are really the NFFA results.

In order to use the FFT, it is necessary to satisfy Eqs. (7.2-19) and (7.2-20). mx 0x and my 0y are the output sampling positions that are modified to do so. It is possible to consider other sampling strategies in which mx 0x and my 0y are replaced by some other output sampling positions, say, xs and ys. The results obtained above are valid in all such cases.

7.5INVERSE DIFFRACTION AND ITERATIVE OPTIMIZATION

Inverse diffraction involves recovery of uðx; yÞ from Uðx0; y0; zÞ as discussed in Chapter 6. In the case of the Fresnel approximation, Eq. (7.2-9) is a Fourier transform, and its inversion yields

ðð

u0ðx; yÞ ¼ U^ðx0; y0; zÞe j2pðfxxþfyyÞdfxdfy

ð7:5-1Þ

NUMERICAL EXAMPLES

97

where u0ðx; yÞ is given by Eq. (7.2-10), the spatial frequencies fx, fy Eqs. (7.2-12), (7.2-13), and

^ð y zÞ ¼ Uðx y zÞe jk2vz

U x0; 0; 0; 0;

are given by

ð7:5-2Þ

Similarly, the integrals in Eqs. (7.3-11) and (7.4-5) are Fourier transforms. Inverse Fourier transformation of Eq. (7.3-11) yields

ðð
u0ðx; yÞ ¼ U0ðx0; y0; zÞe j2pðfxxþfyyÞdxdy	ð7:5-3Þ
where
U0ðx0; y0; zÞ ¼ Uðx0; y0; zÞe jkr0	ð7:5-4Þ
The inverse Fourier transformation of Eq. (7.4-5) similarly yields
ðð
u00ðx; yÞ ¼ U0ðx0; y0; zÞe j2pðfxxþfyyÞdxdy	ð7:5-5Þ

Iterative optimization techniques are often used to design optical devices such as DOE’s [Lee, 1970, 1975, 1979; Zhuang and Ersoy, 1995]. These techniques are discussed in Chapters 15 and 16. For this purpose, many iterations are carried out between the input and output planes. Equations (7.3-11) and (7.5-3) as well as Eqs. (7.4-5) and (7.5-5) are Fourier transform pairs. Hence, iterative optimization can be carried out exactly as before with these equations.

7.6NUMERICAL EXAMPLES

Below the 2-D case is considered first with y variable neglected without loss of generality. In the initial computer simulations, the relevant parameters are chosen as follows:

l ¼ 0:6328 m; x ¼ 1 mm; N ¼ 256; x=l ¼ 7

The minimum z distance for the approximations to be valid as a function of offset x0 at the output plane can be computed by specifying kðr rðapproximateÞÞ to be less than 1 radian. The result with input offset x equal to 1 mm is shown in Figure 7.1. It is observed that for very small ðx0 xÞ, the two approximations behave almost the same, but this quickly changes as ðx0 xÞ increases.

Figure 7.2 shows the ratio z=x0 as a function of the output plane offset x0. This ratio approaches 2 with the new approximation whereas it approaches 36 with the Fresnel approximation as the offset x0 increases. We note that with the FFT, the maximum output offset equals z=ð x=lÞ, x being the input plane sampling