A New Approach to Telescope Fourier Optics

Roy Danchick

Independent Researcher, Los Angeles, CA, USA.

DOI: 10.4236/oalib.1105856                

Abstract

In this paper we take a new approach to solving the Ray-leigh-Sommerfeld-Smythe equation for a telescope’s optical impulse re-sponse to a monochromatic point source spherical wave radiating outward from the origin in the object plane. Two cases are covered: 1) the point source at the origin in the telescope’s near-field object plane; 2) the point source at the origin in the far-field object plane as is the case with satellite infrared sensors when the distance between the telescope and the point source is very much greater than the telescope’s circular aperture diameter. With only the assumption of a thin circular aperture we 1) derive a general solution that erases the distinction among the three classically defined zones of the optical axis: Near, Fresnel, and Fraunhofer; 2) reduce the computational complexity down from two-dimensional to one-dimensional Fourier transform integrals and; 3) identify Filon quadrature as the numerical method of choice for accurately and efficiently approximating the values of these integrals and; 4) provide a computational example.

Keywords

Fourier Optics, Fourier Transforms, Convolution, Bessel Functions, Filon Quadrature

1. Introduction

The Rayleigh-Sommerfeld-Smythe equation is derived from the first principles of Maxwell’s renowned electro-magnetic light wave equations. In this paper we follow the point-source spherical monochromatic wave formulation of the equations in [1] and [2]. In [1] George derives the equation of electric field propagation from the telescope’s circular aperture into its interior. In [2] Dubois derives the simple form of the point-source monochromatic spherical wave equation for the electric field on a telescope’s thin circular aperture.

George makes the traditional distinction between the electric fields in the Near, Fresnel and Fraunhofer zones along the telescope’s optical axis. Then he applies a two-dimensional Fourier analysis approach to approximate the optical impulse response function in the Fraunhofer region. This approximation depends on a linearization after dropping a term that is claimed to be negligible.

In contrast, assuming only a thin circular aperture, we make no distinction between the Near, Fresnel and Fraunhofer zones, keep all terms, do not linearize, and derive reductions from two-dimensional down to one-dimensional Fourier transforms. We cover the two essential cases: 1) the radiating point source delta function at the origin in the near field when the distance between observing sensor and observed object is tens to hundreds of times the telescope diameter; 2) where this distance is very much greater than the telescope’s aperture diameter and is on the order of several tens of millions of meters as is the case with satellite-borne infrared sensor systems.

In both cases the solutions turn out to be the convolution of two circularly symmetric functions of two spatial planar variables. The first of these functions accounts for the propagation of the electric field over the telescope’s aperture into the telescope’s interior while the second is a consequence of the monochromatic spherical wave formulation. The two-dimensional Fourier transform of this first function reduces to a Fourier transform of a function of one spatial variable, the Euclidian distance between the telescope’s interior point and the aperture center, evaluated at the inverse of the wavelength. We derive the one-dimensional Fourier transform of the second function in the near-field first case. In the far-field second case it is the well-known “Sombrero” function.

The rest of the paper is organized as follows: Section 2 details definition and notation. Section 3 is devoted to our derivations. Section 4 provides a computational example. Section 5 draws the conclusions from our derivations and numerical results Section 6 points to possible future research directions. Section 7 supplies our references.

2. Definitions and Notation

Z = Distance along the sensor optical z-axis from the sensor to the point source at the origin of the object plane

${\lambda }_{tel}$ = Operating wave length (m)

$k_{tel}=2\pi /{\lambda }_{tel}$ (m⁻¹) is the optical wave number

$D_{tel}$ = Circular aperture diameter (m)

$f_{tel}$ = Focal length (m)

$\left(x,y,z\right),z>0$ is the coordinate vector of an arbitrary point in the right half space of the telescope’s interior with positive z measured along the optical axis into the telescope from the center of the aperture and the planar x and y coordinates measured orthogonal to the optical z-axis.

The linear transformation between coordinates in the far field object plane, $\left(x_o,y_o,z_o\right)$ and the telescope coordinates $\left(x,y,z\right)$ is given by

$\begin{array}{l}{x}_o=-x,\\ y_o=y,\\ z_o=-z.\end{array}$ (1)

3. Derivations

3.1. The General Form of the Fourier Transform of the Optical Impulse Response Function

Let $E_x\left(x,y,0\right),x^2+y^2\le D_{tel}^2/4$ be the x-component of the electric field on the telescope’s thin circular aperture. Then, following [1], the x-component of the electric field at $\left(x,y,z\right),z>0$ in the telescope interior, propagates from the telescope’s aperture into the telescope interior according to:

$\begin{array}{l}{E}_x\left(x,y,z\right)=\left[iz/\left(2\pi \right)\right]{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)}}\left[\left(i{k}_{tel}+1/R\right)/R^2\right]E_x\left(x^{\prime },y^{\prime },0\right)\text{d}{x}^{\prime }\text{d}{y}^{\prime },\\ R=\sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+z^2}\text{.}\end{array}$ (2)

We first note, treating both Z and z as a parameters, that the integral on the R.H.S of Equation (2) above is a two-dimensional convolution, $P\ast Q$, where

$\begin{array}{l}P\left(x,y;z\right)=\exp \left(-i{k}_{tel}R\right)\left[\left(i{k}_{tel}+1/R^2\right)/R^2\right],\\ R=\sqrt{x^2+y^2+z^2},\\ Q\left(x,y;Z\right)=E_x\left(x,y,0\right)=\exp \left(-i{k}_{tel}S\right)/S,\\ S=\sqrt{x^2+y^2+Z^2}.\end{array}$ (3)

The Fourier transform of Q is, formally:

$\mathbb{F}\left(Q\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp \left[-2\pi i\left(f_{x}x+f_{y}y\right)\right]\exp \left(-i{k}_{tel}S\right)\text{d}x\text{d}y}}/S.$ (4)

If we transform to polar coordinates:

$\begin{array}{l}x=r\cos \theta ,\\ y=r\sin \theta ,\\ r=\sqrt{x^2+y^2},\\ 0\le \theta \le 2\pi ,\\ \text{d}x\text{d}y=r\text{d}r\text{d}\theta ,\\ r\text{d}r=S\text{d}S,\end{array}$ (5)

$\begin{array}{l}\mathbb{F}\left(Q\right)\left(\mu \cos \phi ,\mu \sin \phi \right)\\ ={\displaystyle \underset{x^2+y^2\le D_{tel}^2/4}{\iint }\exp \left[-2\pi i\left(x{f}_x+y{f}_y\right)\right]\left[\exp \left(-i{k}_{tel}S\right)/S\right]\text{d}x\text{d}y}\\ ={\displaystyle \underset{0}{\overset{D_{tel}/2}{\int }}\exp \left(-i{k}_{tel}S\right)r\text{d}r}/S{\displaystyle \underset{0}{\overset{2\pi }{\int }}\exp \left[-2\pi i\mu r\left(\cos \phi \cos \theta +\sin \phi \sin \theta \right)\right]\text{d}\theta }\\ ={\displaystyle \underset{Z}{\overset{\sqrt{Z^2+D_{tel}^2/4}}{\int }}\exp \left(-i{k}_{tel}S\right)\text{d}S}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\exp \left[-2\pi i\mu r\cos \left(\theta -\phi \right)\right]\text{d}\theta }\\ =2\pi {\displaystyle \underset{Z}{\overset{\sqrt{Z^2+D_{tel}^2/4}}{\int }}\exp \left(-i{k}_{tel}S\right)J_0\left(2\pi \mu \sqrt{S^2-Z^2}\right)\text{d}S}\text{,}\end{array}$

$\begin{array}{l}{J}_0\text{is the zeroth order Bessel function of the first kind,}\\ \mu =\sqrt{f_x^2+f_y^2},\\ \phi ={\tan }^{-1}\left(f_y/f_x\right)\text{,}\\ \text{0}\le \phi \le \text{2}\pi \text{.}\end{array}$ (6)

Clearly, the last integral in Equation (6) above exists for all $\mu \in \left[0,\infty \right)$. For practical computational purposes it can be accurately approximated numerically by the method of Filon quadrature [3].

The proof of the existence of $\mathbb{F}\left(P\right)$, which is necessary for both near and far-field case applications is covered below.

3.2. The Fourier Transform of the Far Field Optical Impulse Response Function

We assume $Z\gg D_{tel}$, as is the case for infrared surveillance satellites. We can thus make the highly accurate approximating assumption, which will be validated in the sequel, that the electric field over the circular aperture induced by a 1 watt/steradian point source $\delta$ -function, $\delta \left(x_o,y_o\right)$, monochromatic spherical wave radiating from the origin in the object plane is a constant, $i\exp \left(-i{k}_{tel}Z\right)/Z$. From [1] and [2] we formulate the equation for the electric field impulse response function, $E_x$, to the unit point source, in the telescope’s x-axis direction, onto the telescope aperture and from the aperture into the telescope’s interior point $\left(x,y,z\right),z>0$. An identical treatment holds for the electric field’s y-axis component, $E_y$. The equation for $E_x$ is accurately approximated by:

$\begin{array}{l}{E}_x\left(x,y\right)=\left[i\exp \left(-i{k}_{tel}Z\right)/Z\right]\left[z/\left(2\pi \right)\right]\\ \times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)}}\left[\left(i{k}_{tel}+1/R\right)/R^2\right]\mathrm{circ}\left(2x^{\prime }/D_{tel},2y^{\prime }/D_{tel}\right)\text{d}{x}^{\prime }\text{d}{y}^{\prime },\\ R=\sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+z^2},\\ \mathrm{circ}\left(2x^{\prime }/D_{tel},2y^{\prime }/D_{tel}\right)=\{\begin{array}{l}1\text{if}\sqrt{{\left(x^{\prime }\right)}^2+{\left(y^{\prime }\right)}^2}\le D_{tel}/2\\ 0\text{otherwise}\end{array}\end{array}$ (7)

Again we note, this time treating z as a parameter, that the integral on the R.H.S of Equation (7) above is a two-dimensional convolution, $P\ast Q$ where P is given by Equation (3) and for the far-field case, assuming constant aperture illumination:

$Q\left(x,y\right)=circ\left(2x/D_{tel},2y/D_{tel}\right)$ (8)

If the Fourier transform, $\mathbb{F}\left(P\right)\left(f_x,f_y\right)$, exists we can write the Fourier transform of the convolution, $\mathbb{F}\left(P\ast Q\right)$, as

$\begin{array}{l}\mathbb{F}\left(P\ast Q\right)\left(f_x,f_y\right)=\mathbb{F}\left(P\right)\left(f_x,f_y\right)\mathbb{F}\left(Q\right)\left(f_x,f_y\right),\\ \mathbb{F}\left(Q\right)\left(\mu \cos \phi ,\mu \sin \phi \right))=\left(\pi D_{tel}^2/4\right)J_1\left(\pi D_{tel}\mu \right)/\left(\pi D_{tel}\mu \right),\\ J_1\text{istheBesselfunctionofthefirstkindoforder1,}\end{array}$

$\begin{array}{l}\mu =\sqrt{f_x^2+f_y^2},\\ \phi ={\tan }^{-1}\left(f_y/f_x\right),\\ 0\le \phi \le 2\pi .\end{array}$ (9)

Turning our attention to the existence of $\mathbb{F}\left(P\right)\left(f_x,f_y\right)$, noting the circular symmetry of P as a function of the spatial planar coordinates $\left(x,y\right)$, and once more making the transformation of variables from Cartesian rectangular to polar coordinates, $\mathbb{F}\left(P\right)\left(f_x,f_y\right)$, written formally, is:

$\begin{array}{l}\mathbb{F}\left(P\right)\left(f_x,f_y\right)\\ ={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp \left[-2\pi i\left(x{f}_x+y{f}_y\right)\right]\left[\exp \left(-i{k}_{tel}R\right)\left(i{k}_{tel}+1/R\right)/R^2\right]\text{d}x\text{d}y}}\\ ={\displaystyle \underset{0}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)\left[\left(i{k}_{tel}+1/R\right)/R^2\right]r\text{d}r}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\exp \left[-2\pi i\mu r\left(\cos \phi \cos \theta +\sin \phi \sin \theta \right)\right]\text{d}\theta }\\ ={\displaystyle \underset{0}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)\left[\left(i{k}_{tel}+1/R\right)/R^2\right]r\text{d}r}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\exp \left[-2\pi i\mu r\cos \left(\theta -\phi \right)\right]\text{d}\theta }\\ =2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_0\left(2\pi \mu r\right)\left[\left(i{k}_{tel}+1/R\right)/R^2\right]r\text{d}r},\end{array}$

$\begin{array}{l}{J}_0\text{isthezerothorderBesselfunctionofthefirstkind,}\\ R=\sqrt{x^2+y^2+z^2},\\ r=\sqrt{x^2+y^2},\\ \theta ={\tan }^{-1}\left(y/x\right),\\ 0\le \theta \le 2\pi .\end{array}$ (10)

Thus, we must determine the set of all $\mu$ for which the last integral in Equation (10) above exists for a given $z>0$. We show in the proof below that the integral exists, for a given $z>0$, for all $\mu \in \left[0,\infty \right)$.

Proof

We start by again making a change of variables

$\begin{array}{l}r=\sqrt{R^2-z^2},\\ r\text{d}r=R\text{d}R.\end{array}$ (11)

Let

$\begin{array}{l}u\left(R\right)=J_0\left(2\pi \mu \sqrt{R^2-z^2}\right),\\ \text{d}u=-2\pi \mu J_1\left(2\pi \mu \sqrt{R^2-z^2}\right)R\text{d}R/\sqrt{R^2-z^2},\\ v\left(R\right)=-\exp \left(-i{k}_{tel}R\right)/R,\\ \text{d}v=\exp \left(-i{k}_{tel}R\right)\left(i{k}_{tel}+1/R\right)\text{d}R/R.\end{array}$ (12)

Then, integrating by parts:

$\begin{array}{l}{\displaystyle \underset{0}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_0\left(2\pi \mu r\right)\left[\left(i{k}_{tel}+1/R\right)/R^2\right]r\text{d}r}\\ ={\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_0\left(2\pi \mu \sqrt{R^2-z^2}\right)\left[\left(i{k}_{tel}+1/R\right)/R\right]\text{d}R}\end{array}$

$\begin{array}{l}={\displaystyle \underset{z}{\overset{\infty }{\int }}u\left(R\right)\text{d}v\left(R\right)}={\left[-J_0\left(2\pi \mu \sqrt{R^2-z^2}\right)\exp \left(-i{k}_{tel}R\right)/R\right]}_z^{\infty }\\ -2\pi \mu {\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_1\left(2\pi \mu \sqrt{R^2-z^2}\right)\text{d}R/\sqrt{R^2-z^2}}\\ =\exp \left(-i{k}_{tel}z\right)/z-2\pi \mu {\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_1\left(2\pi \mu \sqrt{R^2-z^2}\right)\text{d}R}/\sqrt{R^2-z^2}.\end{array}$ (13)

We make a second change of variables:

$w^2=\sqrt{R^2-z^2},$ (14)

Then

$\begin{array}{l}2\pi \mu {\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_1\left(2\pi \mu \sqrt{R^2-z^2}\right)\text{d}R}/\sqrt{R^2-z^2}\\ =4\pi \mu {\displaystyle \underset{0}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}z\sqrt{w^4+z^2}\right)J_1\left(2\pi \mu z{w}^2\right)w\text{d}w}/\sqrt{w^4+z^2}.\end{array}$ (15)

From the asymptotic properties of $J_1$ there exist positive numbers, C and T such that for all for all $t\ge T$

$\begin{array}{l}{J}_1\left(t\right)=\sqrt{2/\left(\pi t\right)}\cos \left(t-3\pi /4\right)+c_t/t,\\ \left|c_t\right|\le C.\end{array}$ (16)

Figure 1 below is a graph of the function $t{J}_1\left(t\right)-\sqrt{2t/\pi }\cos \left(t-3\pi /4\right)$. We can see that the choice of $K=0.3$, for example, guarantees that for all $t\in \left(0,\infty \right)$

$\left|J_1\left(t\right)\right|<\sqrt{2/\left(\pi t\right)}+K/t$. (17)

Figure 1. $t{J}_1\left(t\right)-\sqrt{2t/\pi }\cos \left(t-3\pi /4\right)$ vs. t.

From the inequality (17) it follows that there exists a positive number W such that for all $w\ge W$,

$\left|J_1\left(2\pi \mu z{w}^2\right)\right|<\left[1/\left(\pi \sqrt{\mu z}\right)\right]\left[1/w+K/\left(2\sqrt{\mu z}{w}^2\right)\right]$. (18)

Thus

$\begin{array}{l}4\pi \mu \left|{\displaystyle \underset{W}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}z\sqrt{w^4+z^2}\right)J_1\left(2\pi \mu z{w}^2\right)w\text{d}w}/\sqrt{w^4+z^2}\right|\\ \le 4\pi \mu {\displaystyle \underset{W}{\overset{\infty }{\int }}\left|J_1\left(2\pi \mu z{w}^2\right)\right|w\text{d}w}/\sqrt{w^4+z^2}\\ <4\sqrt{\mu /z}{\displaystyle \underset{W}{\overset{\infty }{\int }}\left[1+K/\left(2\sqrt{\mu z}w\right)\right]\text{d}w}/\sqrt{w^4+z^2}<\infty .\end{array}$ (19)

We have now proved that the integrand, $\exp \left(-i{k}_{tel}R\right)J_1\left(2\pi \mu \sqrt{R^2-z^2}\right)/\sqrt{R^2-z^2}$, is absolutely integrable and defines, for each $z>0$, a continuous complex function of $\mu$, $g\left(\mu ;z\right)$, parametrized by z, such that for all $\mu \in \left[0,\infty \right)$ :

$\begin{array}{l}g\left(\mu ;z\right)=\exp \left(-k_{tel}iz\right)/z-2\pi \mu {\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-i{k}_{tel}R\right)J_1\left(2\pi \mu \sqrt{R^2-z^2}\right)\text{d}R}/\sqrt{R^2-z^2},\\ g\left(0;z\right)=\mathrm{e}xp\left(-k_{tel}iz\right)/z.\end{array}$ (20)

If we now define

$f\left(R,\mu ;z\right)=J_0\left(2\pi \mu \sqrt{R^2-z^2}\right)\left(i{k}_{tel}+1/R\right)/R,$ (21)

Equation (11) can be written:

$\begin{array}{l}\mathbb{F}\left(P\right)\left(\mu \cos \phi ,\mu \sin \phi \right)=2\pi {\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;z\right)\text{d}R}=2\pi g\left(\mu ;z\right),\\ 0\le \mu <\infty ,0\le \phi \le 2\pi ,z>0.\end{array}$ (22)

So, $\mathbb{F}\left(P\right)$ exists for all $z>0$ and $\mu \in \left[0,\infty \right)$ and we recognize the integral in Equation (22) above to be the Fourier transform, with respect to R, of the function, $f\left(R,\mu ;z\right)$, evaluated at the spatial frequency $1/{\lambda }_{tel}$. Note that Equation (22) is valid in general on the family of geometric planes parametrized by each $z>0$. So, Equation (22) applies equally in all three classically defined cases: z in the Near, z in the Fresnel, or z in the Fraunhofer zones of the optical axis.

Finally, we can summarize our derivations in the equation for the Fourier transform of the optical impulse response function to a far-field radiating point source, $E_x$, as a function of the polar spatial frequency variables, $\left(\mu ,\phi \right)$ and the parameters Z and z:

$\begin{array}{l}\mathbb{F}\left(E_x\right)\left(\mu \cos \phi ,\mu \sin \phi \right)\\ =\left[iz\exp \left(-2\pi iZ/{\lambda }_{tel}\right)/\left(2\pi Z\right)\right]\mathbb{F}\left(P\right)\left(\mu \cos \phi ,\mu \sin \phi \right)\mathbb{F}\left(Q\right)\left(\mu \cos \phi ,\mu \sin \phi \right)\\ =\left(\pi D_{tel}^2/4\right)\left[iz\exp \left(-2\pi iZ/{\lambda }_{tel}\right)/Z\right]\left[{\displaystyle \underset{z}{\overset{\infty }{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;z\right)\text{d}R}\right]\\ \cdot Somb\left(\pi D_{tel}\mu \right),\\ Somb\left(w\right)=J_1\left(w\right)/w.\end{array}$ (23)

Note that the Equations (14) and (23) above also imply that $\mathbb{F}\left(E_x\right)$ is a continuous function of z, $z\in \left[0,\infty \right)$.

It then follows that, given an object plane intensity-at-the source distribution, $O\left(x_o,y_o\right)$ (watts/steradian), on the far-field object plane $z_o=-Z$, and assuming the existence of its Fourier transform,

$\mathbb{F}\left(O\right)\left(f_x,f_y\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp \left[-2\pi i\left(f_{x}x+f_{y}y\right)\right]O\left(-x,y\right)\text{d}x\text{d}y}}$, (24)

the Fourier transform of the induced electric field in the x-direction, $\mathbb{F}\left(E_{0,x}\right)\left(f_x,f_y\right)$, is given by

$\mathbb{F}\left(E_{0,x}\right)\left(f_x,f_y\right)=\mathbb{F}\left(E_x\right)\mathbb{F}\left(O\right)\left(f_x,f_y\right)$ (25)

and the observed electric field on the geometric plane at z is give by

$E_{0,x}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp \left[2\pi i\left(f_{x}x+f_{y}y\right)\right]\mathbb{F}\left(E_x\right)\mathbb{F}\left(O\right)\left(f_x,f_y\right)\text{d}x\text{d}y}}$. (26)

The issue now becomes one of finding an accurate and efficient numerical method to approximate the Fourier transform integral factor, $\mathbb{F}\left(E_x\right)$, in (25) above for on the geometric focal plane, $z=f_{tel}$, the infinite upper limit replaced by a sufficiently large positive number B so that the magnitude of the integral’s

tail, ${\displaystyle \underset{B}{\overset{\infty }{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;f_{tel}\right)\text{d}R}$ is small relative to

${\displaystyle \underset{f_{tel}}{\overset{B}{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;f_{tel}\right)\text{d}R}$. Equation (22) provides us with two options

for the far-field exploitation of Filon quadrature: 1) apply the Filon method to the integrand $f\left(R,\mu ;f_{tel}\right)$ directly or equivalently 2) apply the Filon method to the integrand that appears in the definition of $g\left(R,\mu ;f_{tel}\right)$. In either case the effectiveness of the numerical method will depend on the behavior of the third and fourth derivatives of these two functions with respect to R. This is the case because the Filon quadrature cosine and sine error terms, $e_{\cos }$ and $e_{\sin }$, can be written for, say, the former case as:

$\begin{array}{c}{e}_{\cos }=C^{\left(3\right)}\left(h\right){\displaystyle \underset{j=1}{\overset{n}{\sum }}{f}^{\left(3\right)}}\left(f_{tel}+\left(2j-1\right)h\right)\sin \left(\left(2j-1\right)j\vartheta \right)\\ -C^{\left(4\right)}\left(h\right)f^{\left(4\right)}\left(\xi \right)+O\left(k_{tel}{h}^7\right),\end{array}$ (27)

$\begin{array}{c}{e}_{\sin }=C^{\left(3\right)}\left(h\right){\displaystyle \underset{j=1}{\overset{n}{\sum }}{f}^{\left(3\right)}}\left(f_{tel}+\left(2j-1\right)h\right)\cos \left(\left(2j-1\right)j\vartheta \right)\\ -C^{\left(4\right)}\left(h\right)f^{\left(4\right)}\left(\xi \right)+O\left(k_{tel}{h}^7\right),\end{array}$ (28)

$C^{\left(3\right)}\left(h\right)=\left(2/45\right)k_{tel}{h}^4,$

$C^{\left(4\right)}\left(h\right)=\left(1/180\right)\left(B-f_{tel}\right)h^4,$

h = the Filon quadrature step-size

$n=\left(B-f_{tel}\right)/h$,

$\vartheta =k_{tel}h$,

$f^{\left(3\right)}$ is the third derivative of f with respect to R,

$f^{\left(4\right)}\left(\xi \right)$ is the fourth derivative of f at the intermediate point $\xi$,

$\xi \in \left(f_{tel},B\right)$.

Equations (27) and (28) above allow us to make two key observations. First, since ${\lambda }_{tel}$ is small, we can choose h to be a large integer multiple of ${\lambda }_{tel}$, say, in our monochromatic case for example, $h={10}^5{\lambda }_{tel}=0.2875\text{m}$. This implies that $\vartheta$ is an integral multiple of 2π and consequently that the first summation term on the R.H.S of (27) for the cosine error term above vanishes identically. We can use a similar trick to make the sine error summation term of (28) to vanish by choosing h to be a large odd integral multiple of ${\lambda }_{tel}/4$. Second, we can choose B to be a large integer multiple of $f_{tel}$, which we can, in turn, choose to be a large integral multiple of ${\lambda }_{tel}$, thus guaranteeing that n is indeed an integer. The decision on which option, (i) or (ii) to choose follows from comparing and contrasting the graphs of the fourth derivatives: ${\text{d}}^{4}f/\text{d}{R}^4$ and ${\text{d}}^4{J}_1\left(2\pi \mu \sqrt{R^2-f_{tel}^2}\right)/\sqrt{R^2-f_{tel}^2}/\text{d}{R}^4$, $R\in \left[f_{tel},B\right]$. Figure 2 and Figure 3 provide the comparative and contrasting graphs that motivate our choice of option (i).

Note how, for $\mu =10$, the magnitude of ${\text{d}}^{4}f/\text{d}{R}^4$ drops off rapidly as R increases, while the magnitude of ${\text{d}}^4{J}_1\left(2\pi \mu \sqrt{R^2-f_{tel}^2}\right)/\sqrt{R^2-f_{tel}^2}/\text{d}{R}^4$ stays much larger than 1 all the way out to B.

4. A Far-Field Computational Example

The idealized geometric setting for this numerical experiment is a geo-stationary

Figure 2. ${\text{d}}^{4}f/\text{d}{R}^4$ vs. R, $\mu =10$.

Figure 3. ${\text{d}}^4{J}_1\left(2\pi \mu \sqrt{R^2-f_{tel}^2}\right)/\sqrt{R^2-f_{tel}^2}/\text{d}{R}^4$ vs. R, $\mu =10$.

infrared surveillance satellite stationed above an ellipsoidal earth at zero degrees latitude and longitude. The telescope’s unit boresight vector is pointed from the satellite to the earth center. The 1 watt/steradian delta function point source is on the equator directly below the spacecraft at zero degrees geodetic latitude and zero degrees longitude. We make the simplifying assumption that there is no atmospheric absorption of the light emitted from the point source so that the source intensity is not attenuated.

4.1. Parameter Values

The numerical values of the parameters are:

$\begin{array}{l}Z=35862945\left(\text{m}\right)\\ {\lambda }_{tel}=2.875\times {10}^{-6}\left(\text{m}\right)\\ k_{tel}=2.1854558\times {10}^6\\ D_{tel}=0.5\left(\text{m}\right)\\ f_{tel}=5.75(m)\end{array}$

4.2. The Electric Field on the Geometric Focal Plane

We focus on the electric field, $E_x$, on the geometric focal plane at $z=f_{tel}$ and, for computational convenience, without loss of generality, have chosen to make the ratio of the focal length to the wavelength an integer,

$f_{tel}/{\lambda }_{tel}=M.$ (29)

We first address the assumption of a uniform electric field on the aperture by comparing the MATLAB computation of the value of the field at the aperture center:

$i\exp \left(-i{k}_{tel}Z\right)/Z=1.561149216428995\text{e}-008-2.310400616893207\text{e}-00\text{8}i$

with the MATLAB computation of the field on the circumference of the aperture:

$\begin{array}{l}i\exp \left(-i{k}_{tel}\sqrt{Z^2+D_{tel}^2/4}\right)/\sqrt{Z^2+D_{tel}^2/4}\\ =1.561149216428995\text{e}-008-2.310400616893207\text{e}-008i\end{array}$

Note that the numerical values of the electric field at aperture center and the field on the circumference agree to 15 decimal places in both their real and imaginary parts. So, we can conclude that the assumption of uniform aperture illumination is indeed valid in this far-field case at hand.

Next we specialize Equation (23) to $z=f_{tel}$, on the geometric focal plane:

Apropos of our analysis of the Filon quadrature fourth order error term we have chosen $B=2001f_{tel}=11505.75\text{m}$ and $h={10}^5{\lambda }_{tel}=0.2875\text{m}$ Table 1 below shows the surprising result that, over the range $\mu \le 1000{\text{m}}^{-1}$ in spatial

frequency the integral, ${\displaystyle \underset{f_{tel}}{\overset{B}{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;f_{tel}\right)\text{d}R}$, is nearly a constant

$\approx 1/f_{tel}$. We note that the small magnitude of its imaginary part is about the size of error term that we see in Figure 2 for $R>3000\text{m}$.

Assuming that the Filon quadrature results shown above are as accurate as they appear to be and hold up as $\mu \to \infty$ we are left with one of two possible inferences: 1) either the Fourier transform,

$\mathbb{F}\left(P\right)\left(\mu \cos \phi ,\mu \sin \phi \right)=2\pi {\displaystyle \underset{f_{tel}}{\overset{\infty }{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;f_{tel}\right)\text{d}R}$, is spatial frequency

band-limited, i.e., has compact support in $\mu \in \left[0,\infty \right)$ so that $\mathbb{F}\left(P\right)$ vanishes outside some very large circle in the spatial frequency domain or 2) begins to approach and tend to zero as $\mu \to \infty$ only for extremely large values of $\mu$.

5. Conclusion

With the single essential hypotheses of a thin circular aperture illuminated by a monochromatic point source at the origin of the object plane (both near-field and far-field) we have proven the existence of and derived exact near and far-field equations for the Fourier transform of the electric field optical impulse response function on any plane perpendicular to the telescope’s optical axis at distance z from the telescope aperture center. This equation applies in general to Near, Fresnel, and Fraunhofer zones of the optical $z$ -axis. Exploiting the method of Filon quadrature to numerically approximate the Fourier transform integral on the geometrical focal plane at $z=f_{tel}$, the results of our far-field numerical experiment show near constancy of the transform for $\mu \in \left[0,1000\right]$. Thus, if indeed, $\mathbb{F}\left(P\right)=2\pi /f_{tel}$ continues to hold for very large $\mu$ and only begins to tend to 0 for extremely large $\mu$, then the approximation (23) with

Table 1. ${\displaystyle \underset{f_{tel}}{\overset{B}{\int }}\exp \left(-2\pi iR/{\lambda }_{tel}\right)f\left(R,\mu ;f_{tel}\right)\text{d}R}$ vs. $\mu$.

$\mathbb{F}\left(P\right)=2\pi /f_{tel}$ would provide a useful optical systems engineering tool.

6. Directions for Further Research

There are three promising areas for future research.

First, is the investigation of the computational accuracy of the Filon method to evaluate $\mathbb{F}\left(Q\right)\left(\mu \cos \phi ,\mu \sin \phi \right)$ in the near-field case. The Filon error Equation (32), applied to the integrand in the near-field Equation (6), would provide the jumping off point for this study.

Second, is the determination whether $\mathbb{F}\left(P\right)\left(\mu \cos \phi ,\mu \sin \phi \right)$ is spatial frequency band-limited in $\mu \in \left[0,\infty \right)$ or tails off to zero as $\mu \to \infty$. A good beginning would be extension of Table 1 to very large values of $\mu$.

Recognizing that our results clearly apply to refractor telescopes, the third area for research is to investigate how our results can be applied to reflector telescope designs such as the Schmidt-Cassegrain telescopes that are used for satellite infrared surveillance missions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	George, N. (2012) Fourier Optics. Monograph, Hajim School of Engineering and Applied Sciences, University of Rochester, Rochester, 8.
[2]	Dubois, A. (2013) Fourier Optics. Monograph, Institut d’Optique, ParisTech, 14.
[3]	Abromowitz, M.A. and Stegun, I. (1972) Mathematical Handbook for Scientists and Engineers. Dover Publications Inc., New York, 890-891.