## Abstract

The numerical calculation of the Rayleigh–Sommerfeld diffraction integral is investigated. The implementation of a fast-Fourier-transform (FFT) based direct integration
(FFT-DI) method is presented, and Simpson's rule is used to improve the calculation accuracy. The sampling interval, the size of the computation window, and their influence on numerical accuracy and on computational complexity are discussed for the FFT-DI and the FFT-based angular spectrum (FFT-AS) methods. The performance of the FFT-DI method is verified by numerical simulation and compared with that of the FFT-AS method.

© 2006 Optical Society of America

Full Article |

PDF Article
### Equations (41)

Equations on this page are rendered with MathJax. Learn more.

(1)
$$\frac{{\partial}^{2}U}{\partial {x}^{2}}+\frac{{\partial}^{2}U}{\partial {y}^{2}}+\frac{{\partial}^{2}U}{\partial {z}^{2}}+{k}^{2}U=0,$$
(2)
$$A\left(\alpha ,\beta ,z\right)=A\left(\alpha ,\beta ,0\right)G\left(\alpha ,\beta ,z\right),$$
(3)
$$A\left(\alpha ,\beta ,z\right)=F\left\{U\left(x,y,z\right)\right\}={\displaystyle \int \int U\left(x,y,z\right)\mathrm{exp}\left(-j\alpha x-j\beta y\right)\mathrm{d}\mathit{}\mathit{x}\mathrm{}\mathrm{d}\mathit{}y},$$
(4)
$$G\left(\alpha ,\beta ,z\right)=\mathrm{exp}\left(j\sqrt{{k}^{2}-{\alpha}^{2}-{\beta}^{2}}z\right)$$
(5)
$$U\left(x,y,z\right)={F}^{-1}\left\{A\left(\alpha ,\beta ,z\right)\right\}$$
(6)
$$=\frac{1}{4{\pi}^{2}}\text{\hspace{0.17em}}{\displaystyle \int \int A\left(\alpha ,\beta ,0\right)\mathrm{exp}\left(j\alpha x+j\beta y+j\sqrt{{k}^{2}-{\alpha}^{2}-{\beta}^{2}}z\right)\mathrm{d}\alpha \mathrm{d\beta}\mathrm{.}}$$
(7)
$$Q=\mathrm{I}\mathrm{F}\mathrm{F}\mathrm{T}2\left\{\mathrm{F}\mathrm{F}\mathrm{T}2\left\{U\left({x}_{m},{y}_{n},0\right)\right\}\text{\hspace{0.17em}}\text{\u22c5 \xd7 \hspace{0.17em}}G\left({\alpha}_{m},{\beta}_{n},z\right)\right\},$$
(8)
$$g\left(x,y,z\right)=\frac{1}{4{\pi}^{2}}{\displaystyle \text{\hspace{0.17em}}\int \int G\left(\alpha ,\beta ,z\right)\mathrm{exp}\left(j\alpha x+j\beta y\right)\mathrm{d}\alpha \mathrm{d}\beta}$$
(9)
$$=\frac{1}{2\pi}\text{\hspace{0.17em}}\frac{\mathrm{exp}\left(jkr\right)}{r}\text{\hspace{0.17em}}\frac{z}{r}\text{\hspace{0.17em}}(\frac{1}{r}-jk),$$
(10)
$$U\left(x,y,z\right)={\displaystyle {\iint}_{A}U\left(\varsigma ,\eta ,0\right)g\left(x-\varsigma ,y-\eta ,z\right)\mathrm{d}\varsigma \mathrm{d}\eta}$$
(11)
$$={\displaystyle {\iint}_{A}U\left(\varsigma ,\eta ,0\right)\text{\hspace{0.17em}}\frac{\mathrm{exp}\left(jkr\right)}{2\pi r}\text{\hspace{0.17em}}\frac{z}{r}(\frac{1}{r}-jk)\mathrm{d\varsigma d}\eta ,}$$
(12)
$$U\left({x}_{m},{y}_{n},z\right)={\displaystyle \sum _{i=1}^{N}{\displaystyle \sum _{j=1}^{N}U\left({\varsigma}_{i},{\eta}_{j},0\right)g\left({x}_{m}-{\varsigma}_{i},{y}_{n}-{\eta}_{j},z\right)}\Delta \varsigma \Delta \eta ,}$$
(13)
$$S=\mathrm{I}\mathrm{F}\mathrm{F}\mathrm{T}2\left[\mathrm{F}\mathrm{F}\mathrm{T}2\left(U\right)\text{\hspace{0.17em}}\text{\u22c5 \xd7 \hspace{0.17em}}\mathrm{F}\mathrm{F}\mathrm{T}2\left(H\right)\right]\Delta \varsigma \Delta \eta ,$$
(14)
$$U={\left[\begin{array}{cc}{U}_{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]}_{(2N-1)\times (2N-1)}=\left[\begin{array}{ccccc}U\left({\varsigma}_{1},{\eta}_{1},0\right)& \cdots & U\left({\varsigma}_{1},{\eta}_{N},0\right)& \text{|}& \\ \vdots \text{}& \ddots & \vdots \text{}& \text{|}& {\text{0}}_{N\times (N-1)}\\ U\left({\varsigma}_{N},{\eta}_{1},0\right)& \cdots & U\left({\varsigma}_{N},{\eta}_{N},0\right)& \text{|}& \\ -& -& -& \text{|}& -\\ & {\text{0}}_{(N-1)\times N}& & \text{|}& {\text{0}}_{(N-1)\times (N-1)}\end{array}\right],$$
(15)
$$H={\left[\begin{array}{ccc}g\left({X}_{1},{Y}_{1},z\right)& \cdots & g\left({X}_{1},{Y}_{2N-1},z\right)\\ \vdots & \ddots & \vdots \\ g\left({X}_{2N-1},{Y}_{1},z\right)& \cdots & g\left({X}_{2N-1},{Y}_{2N-1},z\right)\end{array}\right]}_{(2N-1)\times (2N-1)},$$
(16)
$${X}_{j}=\{\begin{array}{cc}{x}_{1}-{\varsigma}_{N+1-j}& j=1,\dots \text{\hspace{0.17em}},N-1\\ {x}_{j-N+1}-{\varsigma}_{1}& j=N,\dots \text{\hspace{0.17em}},2N-1\end{array},$$
(17)
$${Y}_{j}=\{\begin{array}{cc}{y}_{1}-{\eta}_{N+1-j}& j=1,\dots ,N-1\\ {y}_{j-N+1}-{\eta}_{1}& j=N,\dots ,2N-1\end{array}\mathrm{.}$$
(18)
$$U\left({x}_{m},{y}_{n},z\right)={S}_{m+N,n+N}.$$
(19)
$$E=O\left(\Delta {\varsigma}^{2}\right)+O\left(\Delta {\eta}^{2}\right).$$
(20)
$$U={\left[\begin{array}{cc}W\text{\hspace{0.17em}\u22c5 \xd7}{\text{\hspace{0.17em}}U}_{0}& \mathbf{O}\\ \mathbf{O}& \mathbf{O}\end{array}\right]}_{(2N-1)\times (2N-1)},$$
(22)
$$B=1/3[1\text{\hspace{0.17em}}4\text{\hspace{0.17em}}2\text{\hspace{0.17em}}4\text{\hspace{0.17em}}2\text{\hspace{0.17em}}\dots \text{\hspace{0.17em}}2\text{\hspace{0.17em}}4\text{\hspace{0.17em}}1]$$
(23)
$$E=O\left(\Delta {\varsigma}^{4}\right)+O\left(\Delta {\eta}^{4}\right),$$
(24)
$$k\sqrt{{\left(\rho +\mathrm{\Delta \rho}\right)}^{2}+{z}^{2}}-k\sqrt{{\rho}^{2}+{z}^{2}}=2\pi ,$$
(25)
$$\Delta \rho =\sqrt{{\lambda}^{2}+{\rho}^{2}+2\lambda \sqrt{{\rho}^{2}+{z}^{2}}}-\rho .$$
(26)
$$Q\left({x}_{m},{y}_{n},z\right)={\displaystyle \sum _{i=-\infty}^{\infty}{\displaystyle \sum _{j=-\infty}^{\infty}U\left({x}_{m}+iX,{y}_{n}+jY,z\right)}},$$
(27)
$$\left|g\left(x,y,z\right)\right|<\epsilon \left|g\left(0,0,z\right)\right|,$$
(28)
$$\left|\frac{1}{{r}^{2}}\text{\hspace{0.17em}}(\frac{1}{r}-jk)\left|<\epsilon \right|\frac{1}{{z}^{2}}\text{\hspace{0.17em}}(\frac{1}{z}-jk)\right|,$$
(29)
$$\frac{1}{{r}^{2}}\text{\hspace{0.17em}}<\text{\hspace{0.17em}}\epsilon \text{\hspace{0.17em}}\frac{1}{{z}^{2}},$$
(30)
$$\rho \text{\hspace{0.17em}}>\text{\hspace{0.17em}}z{\left[\left(1/\epsilon \right)-1\right]}^{1/2}.$$
(31)
$${C}_{1}={C}_{a}+{C}_{b}+{C}_{c}+{C}_{d}=O\left({N}^{2}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}N\right)+O\left({N}^{2}\right)+O\left({N}^{2}\right)+O\left({N}^{2}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}N\right)\mathrm{.}$$
(32)
$${C}_{2}={C}_{a}+{C}_{b}+{C}_{c}+{C}_{d}+{C}_{e}$$
(33)
$$=O\left({N}_{F}^{\text{\hspace{0.17em} \hspace{0.17em} \hspace{0.17em}}2}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}{N}_{F}\right)+O\left({N}_{F}^{\text{\hspace{0.17em} \hspace{0.17em} \hspace{0.17em}}2}\right)+O\left({N}_{F}^{\text{\hspace{0.17em} \hspace{0.17em} \hspace{0.17em}}2}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}{N}_{F}\right)+O\left({N}_{F}^{\text{\hspace{0.17em} \hspace{0.17em} \hspace{0.17em}}2}\right)+O\left({N}_{F}^{\text{\hspace{0.17em} \hspace{0.17em} \hspace{0.17em}}2}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}{N}_{F}\right),$$
(34)
$$4\times O\left({N}^{2}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}N\right)=O\left[{N}_{1}^{\text{\hspace{0.17em} \hspace{0.17em}}2}\left({\mathrm{log}}_{2}\text{\hspace{0.17em}}{N}_{1}-1\right)\right],$$
(35)
$$h\left(x,z\right)=\frac{jkz}{2r}\text{\hspace{0.17em}}{H}_{1}^{\text{\hspace{0.17em} \hspace{0.17em}}\left(1\right)}\left(kr\right),$$
(36)
$$U\left(x,z\right)={\displaystyle {\int}_{A}U\left(\varsigma ,0\right)h\left(x-\varsigma ,z\right)\mathrm{d}\varsigma}={\displaystyle {\int}_{A}U\left(\varsigma ,0\right)\text{\hspace{0.17em}}\frac{jkz}{2r}\text{\hspace{0.17em}}{H}_{1}^{\text{\hspace{0.17em} \hspace{0.17em}}\left(1\right)}\left(kr\right)\mathrm{d}\varsigma},$$
(37)
$$S=\mathrm{I}\mathrm{F}\mathrm{F}\mathrm{T}\left[\mathrm{F}\mathrm{F}\mathrm{T}\left(U\right)\text{\hspace{0.17em}}\text{\u22c5 \xd7 \hspace{0.17em}}\mathrm{F}\mathrm{F}\mathrm{T}\left(H\right)\right]\Delta \varsigma ,$$
(38)
$$U={\left[U\left({\varsigma}_{1},0\right)\text{\hspace{0.17em}}\dots \text{\hspace{0.17em}}U\left({\varsigma}_{N},0\right)\text{\hspace{0.17em}}0\text{\hspace{0.17em}}\dots \text{\hspace{0.17em}}0\right]}_{2N-1},$$
(39)
$$H=\frac{jkz}{2}{\left[\frac{{H}_{1}^{\text{\hspace{0.17em} \hspace{0.17em}}\left(1\right)}\left(k{r}_{1}\right)}{{r}_{1}}\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}\frac{{H}_{1}^{\text{\hspace{0.17em} \hspace{0.17em}}\left(1\right)}\left(k{r}_{2N-1}\right)}{{r}_{2N-1}}\right]}_{2N-1},$$
(40)
$${r}_{j}=\{\begin{array}{cc}\sqrt{{\left({x}_{1}-{\varsigma}_{N+1-j}\right)}^{2}+{z}^{2}},& j=1,\dots ,N-1\\ \sqrt{{\left({x}_{j-N+1}-{\varsigma}_{1}\right)}^{2}+{z}^{2}},& j=N,\dots ,2N-1\end{array},$$
(41)
$$U\left(z\right)={U}_{0}z\left[\frac{\mathrm{exp}\left(jkz\right)}{z}-\frac{\mathrm{exp}\left(jk\sqrt{{z}^{2}+{a}^{2}}\right)}{\sqrt{{z}^{2}+{a}^{2}}}\right],$$