Fast Hankel transform and its application for studying the propagation of cylindrical electromagnetic fields

doi:10.1364/OE.10.000521

« Show journal navigation

Fast Hankel transform and its application for studying the propagation of cylindrical electromagnetic fields

D. Zhang, X. Yuan, N. Ngo, and P. Shum »View Author Affiliations

Optics Express, Vol. 10, Issue 12, pp. 521-525 (2002)
http://dx.doi.org/10.1364/OE.10.000521

View Full Text Article

Acrobat PDF (117 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Theory
3. Numerical experiments
4. Application for VPWS
5. Summary
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (5)
Cited By
Figures (4)
Metrics

Abstract

We present a fast Hankel transform (FHTn) method for direct numerical evaluation of electromagnetic (EM) field propagation through an axially symmetric system. Comparing with the vector-based plane-wave spectrum (VPWS) method, we present an alternative approach to implement the fast Hankel transform which does not require an additional coordinate transformation for Fourier transform. The proposed FHTn method is an efficient approach for numerical evaluation of an arbitrary integer order of the Hankel transform (HT). As an example to demonstrate the effectiveness of the proposed method, we apply the FHTn technique to the analysis of cylindrical EM field propagation through a diffractive microlens.

[Optical Society of America ]

1. Introduction

It is known that the plane-wave spectrum (PWS) method or alternatively known as an angular spectrum method, is an efficient technique for studying the propagation of electromagnetic wave from one plane to another by determining the Fourier transform (FT) of the incident field and a multiplicative transfer function over a propagation distance [1

J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, New York, 1968), Chap. 2, pp. 11–13.

]. In numerical analysis of beam propagation through an axially symmetric system, the Fourier transform is replaced by the Hankel transform (HT) with separable variables in polar coordinates.

In 1999, Shi and Prather [2

Shouyuan Shi and Dennis W. Prather, “Vector-based plane-wave spectrum method for the propagation of cylindrical electromagnetic fields,”Opt. Lett. 24, 1445 (1999). [CrossRef]

] reported a vector-based PWS (VPWS) method for the study of the propagation of a cylindrical electromagnetic field that involves reduplicate numerical evaluations of arbitrary integer orders of HT other than the only 0th-orders in many cases. In 1992, Magni et al [3

Vittorio Magni, Giulio Cerullo, and Sandro De Silvestri, “High-accuracy fast Hankel transform for optical beam propagation,”J. Opt. Soc. Am. A. 9, 2031–2033 (1992). [CrossRef]

] proposed an efficient fast Hankel transform algorithm (FHATHA) to implement numerical evaluations of the 0th-order HT. However, the FHATHA method has limited applications in the VPWS method, so it has to require an additional step of coordinate transformation when studying the propagation of the cylindrical wave in the VPWS method. A fast Hankel transform of n-th order has been reported by Ferrari et al [4

José A. Ferrari, Daniel Perciante, and Alfredo Dubra, “Fast Hankel transform of n-th order,”J. Opt. Soc. Am. A. 16, 2581–2582 (1999). [CrossRef]

] in 1999. However, when the order of the transform increases, the sampling frequency of the transform function must be changed in Ferrari’s algorithm. A numerical integration is involved in the algorithm. Since the analytical solution of the integration in Ferrari’s algorithm can not be given directly, numerical evaluation approach should be introduced and thus the computational speed is limited. The FHTn algorithm of this paper does not require the computations of derivatives and numerical integrals, and there is no sampling frequency problem.

In this paper, we present a new approach for the implementation of a fast Hankel transform for direct evaluation of an arbitrary integer order of HT using the VPWS method. The proposed FHTn method is an efficient technique for vector analysis of circular symmetric diffractive optical elements (DOEs) by studying the propagation of an EM field through an arbitrary distance, which can be determined by the finite-difference time-domain (FDTD) method [5

Dennis W. Prather and Shouyuan Shi, “Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements,”J. Opt. Soc. Am. A. 16, 1131–1142 (1999). [CrossRef]

] in the near field region.

As an example, the proposed FHTn method is employed here to perform a numerical experiment to demonstrate its effectiveness in the study of the propagation of an EM field emerging from a diffractive microlens to its focal plane.

2. Theory

It is well known that Fourier transform of a circular symmetric function with separable variables can be written as γth-order HT as follows

g_{γ} (ρ) = \int_{0}^{1} f (r) \cdot J_{γ} (C \cdot rρ) \cdot rdr (0 \leq ρ \leq 1)

(1)

where J_γ is the γth-order Bessel function and C is a constant. We use H_γ {f(r)} to represent the γth-order HT of a function f(r), so Eq. (1) can be re-written as

g_{γ} (ρ) = H_{γ} {f (r)}

(2)

The FHTn method begins with the Bessel function property:

r^{γ + 1} J_{γ} (r) = \frac{d}{dr} [r^{γ + 1} \cdot J_{γ + 1} (r)]

(3)

By substituting Eq. (3) into Eq. (1) with an assumption that the integration is implemented numerically, we can simplify the integral and reduce it to the evaluation of the integral over a small sampling interval. Consequently, in each discrete interval ξ_i < r < ξ _i+1, the value of f(r)/r^γ can be determined by f(r_n )/

r_{n}^{γ}

at a sampling point r _n. Thus the integral over each interval is given by

\int_{ξ_{n}}^{ξ_{n + 1}} f (r) \cdot J_{γ} (C \cdot rρ) \cdot r \cdot dr

= \frac{f (r_{n})}{C \cdot r_{n}^{γ} ρ} \cdot [ξ_{n + 1}^{γ + 1} \cdot J_{γ + 1} (Cρ \cdot ξ_{n + 1}) - ξ_{n}^{γ + 1} \cdot J_{γ + 1} (Cρ \cdot ξ_{n})]

(4)

Following the approach of Ref. [3

Vittorio Magni, Giulio Cerullo, and Sandro De Silvestri, “High-accuracy fast Hankel transform for optical beam propagation,”J. Opt. Soc. Am. A. 9, 2031–2033 (1992). [CrossRef]

], we can rewrite the HT integral of Eq. (1) as a discrete cross-correlation by summing the discrete integrals. The discrete cross-correlation can be numerically evaluated by

g_{γ} (ρ_{m}) = \frac{1}{C ρ_{m}} F {F {φ_{n}} \cdot F^{- 1} {j_{(γ + 1) n}}}

(5)

where F {} represents the Fourier transform (FT) and F ^-1{}the inverse FT. In Eq. (5), the functions j _(γ+1)n and φ_n are defined as

j_{(γ + 1) n} = J_{γ + 1} {C r_{0} \exp [α (n + 1 - N)]}

(6)

for n = 0, 1,…,2N-1,and

φ_{n} = {\begin{matrix} k_{0} F_{0} \exp [α (γ + 1) (1 - N)] for n = 0 \\ F_{n} \exp [α (γ + 1) (n + 1 - N)] for n = 1,2, \dots, N - 1 \\ 0 for = N, N + 1, \dots, 2 N - 1 \end{matrix}

(7)

where f(r_N ) ≡ 0, F_n = f(r_n )/

r_{n}^{γ}

- f(r _n+l)/

r_{n + 1}^{γ}

, k ₀ = [2exp(α)+exp(2α)]/{[1+exp(α)]² [1-exp(-2α)]} , and N is the total number of sampling points. It is seen that the FHTn can be reduced to the FHATHA for the case of γ=0 in Eq. (7). In numerical calculation we only need to reserve the values of g_γ (ρ_m ) from m = 0 to N-1 in Eq. (5) because we pattern the value of φ _n to be zero for n = N, N+1,…,2N-1 according to Eq. (7), therefore the values of g_γ (ρ_m ) from N to 2N-1 are insignificant in mathematics. The parameters r ₀ and α can be chosen arbitrarily in theory, however, for better numerical results, we take the same values of r ₀ =[1+exp(α)]exp(-αN)/2 and α = -ln[1-exp(-α)](N-1) as used in Ref [3

Vittorio Magni, Giulio Cerullo, and Sandro De Silvestri, “High-accuracy fast Hankel transform for optical beam propagation,”J. Opt. Soc. Am. A. 9, 2031–2033 (1992). [CrossRef]

3. Numerical experiments

To test the numerical accuracy of the proposed algorithm, we consider two special functions

Fig. 1. (a)First-order Hankel transform g(y) of f(r)=r (0≤r≤1) for a constant C = 20π ; (b) Difference (error) between the analytical transform and the results obtained with the numerical first-order fast hankel transformations.

Fig.2. (a)Second-order Hankel transform g(y) of _{f(r)=r
⁴
(0≤r≤1)} for a constant C = 20π ; (b) Difference (error) between the analytical transform and the results obtained with the numerical second-order fast hankel transformations.

whose analytical solutions can be easily determined. In this numerical experiment, we implemented the calculations on the Hankel transforms H ₁{r} and H ₂{r ⁴} with a constant C = 20π. The numerical calculations were completed within 0.01 second on a 450MHz PC with 128MByte RAM for 1024 sampling points. Figures 1(a) and 2(a) show the plots of both numerical results g ₁(ρ)and g ₂(ρ). The errors between the numerical results and the analytical solutions are plotted in Figures 1(b) and 2(b). The good agreement between the two results show that the FHTn method is highly accurate, and it can thus be used as a direct Hankel transform approach for studying the propagation through an axially symmetric system.

4. Application for VPWS

We now apply the FHTn method to demonstrate its effectiveness in studying the EM field distribution on the focal plane of a diffractive lens. The lens is normally illuminated with an incident plane wave. For axially symmetric structures, the EM components can be represented in terms of a cylindrical harmonic expansion in the azimuthal direction. From Ref. [2

Shouyuan Shi and Dennis W. Prather, “Vector-based plane-wave spectrum method for the propagation of cylindrical electromagnetic fields,”Opt. Lett. 24, 1445 (1999). [CrossRef]

], the expansion coefficients of the EM components in the focal plane z = z ₀ can be determined with those in the near field z = 0 as follows

E_{ρ}^{m} (l \cdot Δ ρ, z_{0}) = \frac{{(- 1)}^{m + 1}}{2} Δ r^{2}

\times \sum_{i = 0}^{M} i {[E_{r}^{m} (i \cdot Δ r, 0) + E_{ϕ}^{m} (i \cdot Δ r, 0)]

\times k_{0}^{2} H_{m + 1} {J_{m + 1} (k' \cdot l k_{0} Δ ρ) \exp (- j k_{z} z_{0})}

+ [E_{r}^{m} (i \cdot Δ r, 0) - E_{ϕ}^{m} (i \cdot Δ r, 0)]

\times k_{0}^{2} H_{m - 1} {J_{m - 1} (k' \cdot l k_{0} Δ ρ) \exp (- j k_{z} z_{0})}}

(8)

E_{ϕ}^{m} (l \cdot Δ ρ, z_{0}) = \frac{{(- 1)}^{m + 1}}{2} Δ r^{2}

\times \sum_{i = 0}^{M} i {[E_{ϕ}^{m} (i \cdot Δ r, 0) + E_{r}^{m} (i \cdot Δ r, 0)]

\times k_{0}^{2} H_{m + 1} {J_{m + 1} (k' \cdot l k_{0} Δ ρ) \exp (- j k_{z} z_{0})}

+ [E_{ϕ}^{m} (i \cdot Δ r, 0) - E_{r}^{m} (i \cdot Δ r, 0)]

\times k_{0}^{2} H_{m - 1} {J_{m - 1} (k' \cdot l k_{0} Δ ρ) \exp (- j k_{z} z_{0})}}

(9)

E_{z}^{m} (l \cdot Δ ρ, z_{0}) = {(- 1)}^{m + 1} \cdot {(Δ r)}^{2}

\times \sum_{i = 0}^{M} i \cdot {E_{z}^{m} (i \cdot Δ r, 0)

\times k_{0}^{2} H_{m} {J_{m} (k' \cdot l k_{0} Δ ρ) \times \exp (j k_{z} \cdot z_{0})}

(10)

where j = √-1 ,

k_{z} = {[k_{0}^{2} - (k' k_{0})]}^{\frac{1}{2}},

, k ₀ = 2π/λ, 0 ≤ k’ ≤ 1, Δr and Δρ are the sampling lengths in the original and propagated planes, respectively, and λ is the wavelength of illumination light. The index i and l represent the values of the i-th sampling point in the near field and the l-th in the focal plane, respectively. M is the total number of sampling points in the near field, and the constant C in the process of HT is equal to Nk ₀Δr.

In this example, since the incident wave is a normal illumination, only the m = 1 mode is needed. Thus, only zero-, first-, and second-order HTs are used in the calculation. For oblique incidence, higher-order HT will be involved in the VPWS calculation, and in this case, the FHTn method can still be used because it is a generic formula for an arbitrary integer order HT.

For convenience of comparison, we used the same eight-level diffractive lens that was studied in Ref. [5

]. The lens parameters include relative permittivity ∊_r = 2.25; wavelength λ = 1.0 μm; diameter = 102.47 μm; f-number = 0.78; and 15 zones. Firstly, we employed our axially symmetric finite-difference time-domain (FDTD) program to obtain the near electric field and then used Eqs. (8) to (10) to calculate the electric field propagated onto the focal plane. Figure 3 shows the near electric field magnitude at the emergent plane of the lens as determined by the axially symmetric FDTD method. Figure 4 shows the numerical result of the electric field distribution at the focal plane of the lens, which is the same as that determined with the EM method presented in Ref. [5

]. For this example, the VPWS method with FHTn took 41s and only 1M byte to propagate, compared with 913s and 897M byte for the conventional EM method. It should be noted that the proposed FHTn method is an efficient alternative to implement the VPWS analysis since no coordinate transformations are required in the calculation.

Fig. 3 Electric field magnitude at the emergent plane of the lens.

Fig. 4 Electric field magnitude in the focal plane of the lens.

5. Summary

We have presented a fast arbitrary integer order FHTn method for direct numerical evaluation of the field propagation through an axially symmetric system. It is an efficient and fast technique for evaluating arbitrary integer order HT. Comparing with Prather’s VPWS method, the proposed method does not require a coordinate transformation for an axially symmetric problem.

References and Links

1.	J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, New York, 1968), Chap. 2, pp. 11–13.
2.	Shouyuan Shi and Dennis W. Prather, “Vector-based plane-wave spectrum method for the propagation of cylindrical electromagnetic fields,”Opt. Lett. 24, 1445 (1999). [CrossRef]
3.	Vittorio Magni, Giulio Cerullo, and Sandro De Silvestri, “High-accuracy fast Hankel transform for optical beam propagation,”J. Opt. Soc. Am. A. 9, 2031–2033 (1992). [CrossRef]
4.	José A. Ferrari, Daniel Perciante, and Alfredo Dubra, “Fast Hankel transform of n-th order,”J. Opt. Soc. Am. A. 16, 2581–2582 (1999). [CrossRef]
5.	Dennis W. Prather and Shouyuan Shi, “Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements,”J. Opt. Soc. Am. A. 16, 1131–1142 (1999). [CrossRef]

OCIS Codes
(050.1970) Diffraction and gratings : Diffractive optics
(070.2580) Fourier optics and signal processing : Paraxial wave optics
(090.1970) Holography : Diffractive optics

ToC Category:
Research Papers

History
Original Manuscript: May 24, 2002
Revised Manuscript: June 11, 2002
Published: June 17, 2002

Citation
D. Zhang, X. Yuan, N. Ngo, and P. Shum, "Fast Hankel transform and its application for studying the propagation of cylindrical electromagnetic fields," Opt. Express 10, 521-525 (2002)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-10-12-521

Sort: Journal | Reset

References

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper