Boundary diffraction wave integrals for diffraction modeling of external occulters

Eric Cady

doi:10.1364/OE.20.015196

1. Introduction

One of the major goals of the coming decades is to directly image a terrestrial exoplanet. Direct imaging allows spectral information about the planet’s atmosphere to be extracted from the spectrum of reflected light, including the potential presence of biomarkers, gases whose presence and abundance may indicate the presence of life. Detecting planets directly is a difficult task, however, for two reasons: intensity ratio (contrast) and angular resolution. An Earth-twin which orbits a Sun-twin at 1AU emits 10¹⁰ times less flux than the parent star in the visible band; if the system was located at 10 parsecs from Earth, the angular separation of the two objects would be 100mas, which for most proposed space telescopes puts the separation under 4λ/D. Specialized methods are required to remove this flux at these small working angles.

One proposed method of suppressing the starlight is an occulter, a spacecraft with a shaped edge flown in front of the telescope to block the starlight before it arrives at the telescope. The occulter size (tens of meters) and distance (tens of thousands of kilometers) are chosen so the angular extent of the occulter is smaller than some desired angle, on the order of 100 milliarcseconds, so exoplanets in orbit about the star will still be visible. The edge is shaped to control diffraction, with the form chosen to suppress the light to a factor of 10¹⁰ across the telescope aperture and over a wide spectral passband.

A major advantage of this approach, compared to an internal coronagraph, is that it significantly loosens the tolerances on the optical quality and thermal stability of the telescope optics, as well as eliminating the need for a wavefront control system for active correction. Conversely, the occulter itself must be built, deployed, and flown in formation to its own set of tolerances, and demonstrating this capability is an area of active research [1]. These tolerances include limits on permanent errors in the occulter shape due to manufacturing and deployment errors; transient shape deformations from thermal fluctuations and vibration; and position and orientation changes of the occulter relative to the telescope-target axis from formation flying [2, 3].

The scale of an occulter—and the distance it must maintain with respect to its target—precludes optical testing on the ground. All of these effects must be modeled in order to build an error budget and verify that an occulter can satisfy these requirements. The capture of the large dynamic range in the resulting field—10⁵ or greater in amplitude—drives the accuracy of the propagator used for the models, while the modeling of time-varying thermal shape deformations and closed-loop formation flying in broadband light drive the speed of the calculation. Some of the most in-depth simulated observations can take from hours to days to run, and the ability to quickly and accurately model propagation past an occulter is thus a major enabling factor in the characterization of occulter system performance. Section 2 gives an overview of the current techniques; the new technique presented in this work is described in Sections 3 and 3.1, and a suggested implementation scheme, along with a pseudocode overview of the algorithm and some computational results, are given in Section 4.

2. Current occulter modeling techniques

Occulter designs take the general form of a solid central disk surrounded by N_p identical tapering structures (petals) around the edge, giving the whole structure the general appearance of a flower. (An example is shown in Fig. 1.) This shape is a result of their provenance as (0,1)-valued approximations to apodizers; the general method to design an occulter is to determine a smooth apodization profile A(r) which can provide the necessary starlight suppression, and then convert it into a binary shape with a sufficient number of petals so that performance is not undermined. The actual form of A(r) can be determined by optimization [4, 5] or selecting an existing functional form [6, 7].

Fig. 1 An example of a typical occulter with 20 petals, with the associated Cartesian and polar coordinate systems.

Download Full Size | PDF

For an occulter designed in this manner, the region Ω in which the occulter is opaque can be written directly:

Ω = {(r, θ) : 0 \leq r \leq R, θ \in Θ (r)}

where Θ (r) = \cup_{n = 0}^{N_{p} - 1} [\frac{2 π n}{N_{p}} - \frac{π}{N_{p}} A (r), \frac{2 π n}{N_{p}} + \frac{π}{N_{p}} A (r)],

with (r, θ) being polar coordinates in the plane of the occulter. Each of the [...,...] is a set of points which defines the outline of a single petal; ∪ denotes the union of these sets, which sets the boundary for the entire occulter. (We denote this boundary as ∂Ω, as it will be used later.)

For a typical occulter, the occulter-telescope distance z = 5 × 10⁷m; the wavelength λ = 5 × 10⁻⁷m; the occulter radius R = 25m; and the maximum excursions in the plane of the telescope |ξ| and |η| ≤ 3m. These values place the occulter well within the Fresnel regime; that is, the region in which it is valid to approximate the exponent in the propagation integral with the first two terms of a power series (the Fresnel approximation). This approximation is suitable when [8].

z^{3} ≫ \frac{π}{4 λ} {[{(x - ξ)}^{2} + {(y - η)}^{2}]}_{max}^{2} .

Further, a standard propagation integral computes the electric field past a finite aperture; the occulter, on the other hand, is a finite obstruction with an effectively infinite aperture. Here, we can turn to Babinet’s theorem [9], which notes that field propagated through an aperture and the field propagated past its complement sum to the field that would have resulted were no obstruction present at all.

If a plane wave with amplitude A and wavelength λ is normally-incident on the occulter, as would be the case for an occulter properly aligned with a target star, the electric field at a distance z past the occulter would be written with the Fresnel approximation and Babinet’s principle as

\begin{array}{l} U (ξ, η) & = & A exp (\frac{2 π i z}{λ}) (1 - \frac{1}{i λ z} \int \int_{Ω} exp {\frac{π i}{λ z} [{(r cos θ - ξ)}^{2} + {(r sin θ - η)}^{2}]} r d r d θ) \\ = & A exp (\frac{2 π i z}{λ}) (1 - \frac{1}{i λ z}) \int \int_{Ω} exp {\frac{π i}{λ z} [{(x - ξ)}^{2} + {(y - η)}^{2}]} d x d y), \end{array}

where ξ and η are Cartesian coordinates at the downstream location. The “1–” term before the integral results from the use of Babinet’s principle. Errors in (for example) occulter shape will modify Ω, but Eq. (4) will hold regardless. Computing the integral for points in the plane of a telescope aperture is the key step in modeling the performance of an occulter.

One could attempt this integral by brute force: discretizing the shape of the occulter over a sufficiently fine grid to capture all of the necessary structure, and propagating this to the telescope aperture with direct integration or a matrix Fourier transform[10]. This approach is not an easy task, as the occulter will be tens of meters across, with tolerances at the level of tens of microns[2], requiring a grid up to 10⁶ × 10⁶ or more to capture the full shape of the profile, not including array padding.

A better approach is to take advantage of the structure of the occulter: in the absence of errors, it has an N_p-fold symmetry in the azimuthal direction, and most errors introduce small perturbations to that. Taking advantage of the symmetry leads to three main approaches to modeling propagation past an external occulter:

Integrals in r: The Bessel function expansion [4] does the integral over θ explicitly, reducing the 2D integral to a series of single integrals in r whose contributions fall off exponentially fast near the optical axis, which makes computation extremely quick in this region. The series can be modified to incorporate some shape errors [11] at the expense of slower convergence. Unfortunately, it is slower at locations far from the optical axis, and requires a good deal of effort to incorporate many types of errors, such as occulter tilt, into the model. Aside from the propagation modeling, the first, radially-independent term in this series is used in optimization approaches [4, 5] to design occulter shapes.
Integrals in θ: The Dubra-Ferrari integral [12] starts one step back from the Fresnel form of Eq. (4), at the Raleigh-Sommerfeld diffraction integral, and evaluates the integral over r to produce a pair of single integrals in angular coordinates, one of which holds for points inside the geometric extent of the starshade, and one which holds outside. This approach can incorporate most errors straightforwardly and spends the same amount of time to determine the field any point in the telescope aperture plane. Because of the manner in which the integrals are segmented, however, the integrands are more difficult to determine correctly for points physically outside the extent of the starshade.
Perturbations to a nominal shape: The slit approach [13] assumes the electric field for an occulter with no errors has been calculated by one of the above methods, and approximates the difference between the shape with and without errors as a series of long, thin boxes whose Fresnel transforms can be calculated quickly and exactly. This method is very fast but necessarily approximate.

In this paper, we propose a fourth method—the boundary-diffraction-wave approach—which seems to out-perform the three extant methods for general purpose use. It maintains the constant-time property of the Dubra-Ferrari integral with a 2–3× improvement in computational speed; see Sec. 4.1 for further discussion. It incorporates all errors on the shape of the occulter natively, and can be written straightforwardly to include position and orientation errors, as well as model waves from off-axis sources, such as an exoplanetary system or the finite diameter of the star. It can also take advantage of the structure of the integral to evaluate the field at multiple wavelengths with little additional computation. The approach of defining the occulter solely in terms of its boundary is a good fit for tolerancing the sorts of errors that are expected for an occulter. Manufacturing errors can be applied by specifically deforming desired edge regions; errors in petal placement can be included by translating and rotating only the boundary points associated with that petal.

3. The boundary diffraction wave formulation

Most diffraction approaches to propagation through a finite aperture find their basis in the Huygens-Fresnel principle, which expresses the electric field at a point following the aperture as the sum over spherical waves emitting from every point within the aperture. However, a different representation was introduced by Maggi and Rubinowicz which splits the field into two parts, a part which propagates according to geometrical optics and a part which incorporates the effect of diffraction from the boundary. This boundary-diffraction-wave (BDW) representation was codified by the work of Miyamoto and Wolf [14, 15], who showed that the boundary integral could be represented for any incident field entirely by a line integral around a vector potential. We will review their analysis briefly here.

Miyamoto and Wolf [14] begin by writing the Kirchhoff formulation (see e.g. [9]) for the diffraction integral from the aperture. Consider a volume of space, bounded by a surface S; the field at any point P within the volume may be expressed as an integral over all surface points Q on S:

U (P) = \frac{1}{4 π} \iint_{S} (U (Q) \nabla_{Q} \frac{exp (i k s)}{s} - \frac{exp (i k s)}{s} \nabla_{Q} U (Q)) \cdot n d S \equiv \frac{1}{4 π} \iint_{S} V (Q, P) \cdot n d S .

Here we let s be the distance between P and Q, n be a vector normal to S at point Q, and ∇_Q is the gradient evaluated at point Q. In the case of the occulter, the volume is the space z ≥ 0, and S stretches over the z = 0 plane, as well as having a component at infinity which is assumed to vanish. (Showing that this term vanishes requires some minor additional assumptions; see Sec. 8.3.2 of [9]).

It can then be shown [14] that for V of this form, there exists a vector potential W, such that V = ∇ × W, and for plane waves incident on an aperture, that associated vector potential is

W (Q, p) = \frac{1}{4 π} A exp (i k p \cdot r^{'}) \frac{exp (i k s)}{s} \frac{\hat{s} \times p}{1 + \hat{s} \cdot p},

U (P) = \iint_{S} \nabla \times W (Q, P) \cdot n d S .

Here r is the vector from the origin to P; Q is a point within the aperture on S, and r′ is the vector from the origin to Q. k = 2π/λ is the wavenumber and λ is the wavelength under consideration. The various s-variables become

s \equiv r^{'} - r,

s = ‖ s ‖,

\hat{s} = \frac{s}{s} .

and we assume the form of the plane wave in the vector direction p to be

U_{0} (P) = A exp (i k p \cdot r),

with amplitude A. A coordinate diagram is shown in Fig. 2.

Fig. 2 A diagram of the points, coordinate systems, and vectors used in this paper. The left grid is in the plane of the occulter, and the right grid in the plane of the telescope aperture.

Download Full Size | PDF

Next, we consider for the moment the case of an opaque screen with an opening over the region Ω, a subregion of S with boundary ∂Ω. We recall here that S is the entire bounding surface of the half-plane z ≥ 0, including the z = 0 plane; Ω is the section of it on which the occulter lies. (We are not considering the occulter case at the moment, but its complement; we will return to the occulter shortly.) In this case, the electric field from Eq. (7) across the aperture is:

U_{ap} (P) = \iint_{Ω} \nabla \times W (Q, P) \cdot n d Ω,

Our approach is to reduce this integral to a line integral, using Stokes’ Theorem. To do this, we consider singularities of the vector potential.

In general, W will have singularities somewhere on S for a given P, as otherwise the field in the half-space past the aperture will be zero. In the specific case of the plane wave, the lone singularity occurs at 1 + ŝ· p = 0, which occurs only when ŝ is parallel to the propagation direction of the plane wave; thus, the singularity will fall in the aperture only for points P which fall in the beam as dictated by geometric optics. (Fig. 3 shows an example of this.) Miyamoto and Wolf [15] show that when Stokes’ Theorem is applied, we get two distinct cases for the field at point P, depending on the singularity location:

\begin{array}{l} U_{ap} (P) & = & A exp (i k p \cdot r) + U^{(B)} (P) for points inside the extent of the aperture \\ = & U^{(B)} (P) otherwise, \end{array}

where U^(B)(P) is the counterclockwise line integral about the edge of Ω:

U^{(B)} (P) = \int_{\partial Ω} W (Q, P) \cdot ℓ d ℓ .

Here, ℓ is a unit vector in the direction tangent to the edge at any point, and dℓ is a differential element of the boundary.

Fig. 3 A diagram showing the projected extent of the aperture given by geometric optics. A plane wave incident on the aperture in plane A with vector direction p produces a shifted aperture in plane B. Two angles which parameterize the direction of the plane wave, ψ₁ and ψ₂, are shown as well.

Download Full Size | PDF

This formulation gives the field following an aperture, but in our case we have the opposite: the occulter is opaque, and the surroundings are open. Here we return to Babinet’s theorem. As we started with the plane wave in Eq. (11), this implies our occulter field U(P) satisfies

U_{0} (P) = U (P) + U_{ap} (P),

or

\begin{array}{l} U (P) & = & - \int_{\partial Ω} W (Q, P) \cdot ℓ d ℓ for points in the geometric shadow of the occulter, \\ = & A e^{i k p \cdot r} - \int_{\partial Ω} W (Q, P) \cdot ℓ d ℓ otherwise . \end{array}

The part of Eq. (16) which depends on geometric optics is termed U^(G)(P):

\begin{array}{l} U^{(G)} (P) & = & 0 for points inside the extent of the aperture \\ = & U_{0} (P) otherwise . \end{array}

This equation is the general form of the occulter field, in the boundary-diffraction-wave formulation.

3.1. Vector potentials for occulter propagation

From this suitably general representation, we can make appropriate substitutions for the specific case of an occulter; these will simplify computation of electric fields later. We define the vector from the origin to P as r = (ξ, η, z); r′, the vector from the origin to Q, is defined as r′ = (x, y, 0). (See again Fig. 2.) We can write the s-variables explicitly as

s \equiv r^{'} - r = (x - ξ, y - η, - z),

s = ‖ s ‖ = {[z^{2} + {(x - ξ)}^{2} + {(y - η)}^{2}]}^{1 / 2},

\hat{s} = \frac{s}{s} .

We consider first the case where the plane wave is normally-incident on the occulter, as starlight would be when the system is correctly aligned. Given this, our vector terms can be written as

p = (0, 0, 1),

p \cdot r^{'} = 0,

\hat{s} \times p = \frac{1}{s} (y - η, - (x - ξ), 0),

1 + \hat{s} \cdot p = 1 - \frac{z}{s},

and the potential in Eq. (6) becomes

exp (i k s) = exp (i k z) exp [i k (s - z)],

W (Q, P) = \frac{1}{4 π} A exp (i k z) \frac{exp [i k (s - z)]}{s - z} \hat{s} \times p .

More generally, an off-axis source at an angle of ψ₁ off-axis and ψ₂ from the x-axis will have

p = (sin ψ_{1} cos ψ_{2}, sin ψ_{1} sin ψ_{2}, cos ψ_{1}),

p \cdot r^{'} = x sin ψ_{1} cos ψ_{2} + y sin ψ_{1} sin ψ_{2} .

Before the expression for the vector potential becomes too complicated, we will define three intermediate variables (f, g, and h) which we can substitute into Eq. (19); these choices will greatly simplify subsequent computations.

f = z sin ψ_{1} + (x - ξ) cos ψ_{2} cos ψ_{1} + (y - η) sin ψ_{2} cos ψ_{1},

g = - (x - ξ) sin ψ_{2} + (y - η) cos ψ_{2},

h = z cos ψ_{1} - (x - ξ) cos ψ_{2} sin ψ_{1} - (y - η) sin ψ_{2} sin ψ_{1} \equiv s \cdot p .

Substituting these variables gives:

s = {[f^{2} + g^{2} + h^{2}]}^{1 / 2},

1 + \hat{s} \cdot p = 1 - \frac{h}{s},

\begin{array}{l} \hat{s} \times p & = & \frac{1}{s} (f sin ψ_{2} + g cos ψ_{2} cos ψ_{1}, \\ - f cos ψ_{2} + g sin ψ_{2} cos ψ_{1}, \\ - g sin ψ_{1}) \end{array}

\begin{array}{l} exp (i k h) & = & exp [i k z cos ψ_{1} - i k (x - ξ) cos ψ_{2} sin ψ_{1} - i k (y - η) sin ψ_{2} sin ψ_{1}] \\ = & exp (- i k p \cdot r^{'}) exp (i k z cos ψ_{1}) exp [i k (ξ cos ψ_{2} + η sin ψ_{2}) sin ψ_{1}], \end{array}

W (Q, P) = \frac{1}{4 π} A exp (i k z cos ψ_{1}) exp [i k (ξ cos ψ_{2} + η sin ψ_{2}) sin ψ_{1}] \frac{exp (i k (s - h))}{s - h} \hat{s} \times p,

a similar relation to Eq. (26). This equation is still exact; for ψ₁ ≪ 1, the usual case for objects in an exoplanetary system, we can note that

\frac{1}{s (s - h)} = \frac{1}{f^{2} + g^{2}} (1 + \frac{1}{{[1 + \frac{f^{2} + g^{2}}{h^{2}}]}^{1 / 2}}) \approx \frac{2}{f^{2} + g^{2}},

s - h \approx \frac{f^{2} + g^{2}}{2 h},

and so the potential becomes

W (Q, P) \approx \frac{1}{2 π} A exp (i k z cos ψ_{1}) exp [i k (ξ cos ψ_{2} + η sin ψ_{2}) sin ψ_{1}] \frac{exp (i k \frac{f^{2} + g^{2}}{2 h})}{f^{2} + g^{2}} v (f, g),

\begin{array}{l} v (f, g) & = & (f sin ψ_{2} + g cos ψ_{2} cos ψ_{1}, \\ - f cos ψ_{2} + g sin ψ_{2} cos ψ_{1}, \\ - g sin ψ_{1}) = s \times p . \end{array}

This form proves to be particularly useful, as (f² + g²) and v · ℓ can be calculated without loss of precision, as neither will be of order z for the small ψ₁ case. A typical occulter might have z = 5×10⁷m, R = 25m, and |ξ| and |η| ≤ 3m, and k ≈ 10⁷m⁻¹; a representative exoplanet target for this occulter might have ψ₁ = 5 × 10⁻⁷ rad. We note that exponential terms of order kz have been separated from smaller terms; these should be calculated independently, as kz ∼ 10¹⁴, and combining terms will lose precision in evaluating the exponent.

4. Efficient implementation

Evaluating the field at each downstream point requires the determination of two parts: the geometric field and the boundary field. We assume, to begin, that we are given a set of N points representing the edge of the occulter which form a simple closed curve, as well as a list of M downstream points and a list of L wavelengths at which the field is to be determined. We also assume that the edge is sufficiently well-sampled that the section of the edge between explicitly-specified points is well-modeled by a linear segment, as this allows us to use the midpoint rule for numerically approximating the integral. Given the slowly-changing shape over the majority of the petal, this is a good assumption; small regions such as petal tips may be specified more finely.

The geometric field is still a plane wave of form Eq. (11) outside the aperture, and is 0 everywhere inside; finding the geometric field thus becomes a problem of checking whether a downstream point falls behind it or not. We do this by treating the edge as a polygon and finding the winding number of the downstream point (ξ₀,η₀); this can be done efficiently with an O(N) routine such as polywind in Numerical Recipes [16], and holds regardless of the complexity of the occulter shape. This approach also speeds up multiband calculations, as the geometric extent of the shadow determined this way is wavelength-independent. Lateral errors in the occulter-telescope position may be included by adding a Δξ and Δη to all points at the telescope plane, and that should be done first. Note that in the case of an off-axis source, the geometric shadow is shifted laterally as well; in this case, the winding number should be calculated around (ξ₀ − zsinψ₁ cosψ₂,η₀ − zsinψ₁ sinψ₂). We then have:

\begin{array}{l} U^{(G)} (ξ_{0}, η_{0}, z) & = & A exp (i k z cos ψ_{1}) exp [i k (ξ_{0} cos ψ_{2} + η_{0} sin ψ_{2}) sin ψ_{1}], \\ for points outside the shadow \\ = & 0 otherwise . \end{array}

In some cases—for example, when the telescope aperture remains close to the center of the occulter and the perturbations to the ideal shape of the occulter are small—the winding number may be able to be determined independently. This determination should be done if possible, as even an efficient algorithm can increase runtime significantly for boundaries containing a large number of points.

The boundary field can be evaluated efficiently with midpoint-rule quadrature running around the edge. First, for each of the N line segments running counterclockwise around the occulter edge, we derive the midpoint of line segment j going from (x_j, y_j, 0) to (x_j₊₁, y_j₊₁, 0):

(x_{j}^{(m)}, y_{j}^{(m)}, 0) = ([x_{j} + x_{j + 1}] / 2, [y_{j} + y_{j + 1}] / 2, 0),

and the vector pointing along the segment

{(x_{j}^{(ℓ)}, y_{j}^{(ℓ)}, 0)}_{j} = (x_{j + 1} - x_{j}, y_{j + 1} - y_{j}, 0) .

with (x_j+1, y_j+1) equal to (x₁, y₁) when j = N. This derivation need only be done once at the beginning, as it holds for all M points and L wavelengths. (We note that, while the midpoint rule has been used for the quadrature, this is certainly not a requirement, and higher order methods could be used; Eq. (39) will still hold. Using a higher-order method may come at the cost of runtime.)

For each point in the downstream field, we calculate f, g, and h for all N segments following Eq. (29):

f_{j} = z sin ψ_{1} + (x_{j}^{(m)} - ξ_{0}) cos ψ_{2} cos ψ_{1} + (y_{j}^{(m)} - η_{0}) sin ψ_{2} cos ψ_{1},

g_{j} = - (x_{j}^{(m)} - ξ_{0}) sin ψ_{2} + (y_{j}^{(m)} - η_{0}) cos ψ_{2},

h_{j} = z cos ψ_{1} - (x_{j}^{(m)} - ξ_{0}) cos ψ_{2} sin ψ_{1} - (y_{j}^{(m)} - η_{0}) sin ψ_{2} sin ψ_{1},

and we then need only two intermediate terms:

T_{1 j} = f_{j}^{2} + g_{j}^{2},

\begin{array}{l} T_{2 j} & = & v (f_{j}, g_{j}) \cdot ℓ d ℓ \\ = & x_{j}^{(ℓ)} (z sin ψ_{1} sin ψ_{2} + (y_{j}^{(m)} - η_{0}) cos ψ_{1}) - y_{j}^{(ℓ)} (z sin ψ_{1} cos ψ_{2} + (x_{j}^{(m)} - ξ_{0}) cos ψ_{1}), \end{array}

to give us:

\begin{array}{l} U^{(B)} (ξ_{0}, η_{0}, z) & = & \frac{1}{2 π} A exp (i k z cos ψ_{1}) exp [i k (ξ_{0} cos ψ_{2} + η_{0} sin ψ_{2}) sin ψ_{1}] \\ \times \sum_{j = 1}^{N} \frac{exp (i k \frac{T_{1 j}}{2 h_{j}})}{T_{1 j}} T_{2 j} . \end{array}

Note that all of the steps prior to Eq. (49) are wavelength-independent; different values of k may be iterated over in the final step. It is advisable to keep terms of order kz (e.g. kzcosψ₁ or kh) separate when evaluating to maintain numerical precision. A typical implementation might take the form shown in Algorithm 1. The for-loops over j in particular are well-suited to vectorization where this is supported.

Algorithm 1:. Pseudocode representation of BDW calculation

View Table

4.1. Computational results

In the standard formulation of the Dubra-Ferrari integral (DF), points outside the geometric shadow are calculated with a different set of integrals which require searches for nearest- and further-neighbor points along the boundary for each point at the telescope aperture. To simplify direct comparison between the Dubra-Ferrari approach and the boundary-diffraction-wave approach for this paper, we choose to select a region which lies entirely in the geometric shadow of the occulter. This selection also fixes the winding number to 1 at all points under consideration, letting us bypass an explicit calculation. This choice is not an unreasonable assumption; the telescope will generally be entirely in the geometric shadow for a properly aligned occulter-telescope system.

Both approaches have similar formulations, with the primary difference coming in the fact that the Dubra-Ferrari integral uses angular coordinates θ centered about (ξ₀, η₀), and so the angular distance between occulter edge points must be rederived for each point at the telescope aperture. This calculation turns out to be one of the major sources of overhead during the DF calculation.

To compare the runtimes of the Dubra-Ferrari and BDW implementations, we consider a closed occulter boundary specified by 164000 individual points, corresponding to the occulter shown in Fig. 1. We select subsets of these points, and run Matlab implementations of both algorithms on this subset. Figure 4 shows the time to implement the propagation as a function of number of boundary points; the BDW algorithm remains 2–3× faster over a wide range of polygon sizes. Figure 5 shows the intensity profile in a region around the telescope aperture, plotted on a base-10 logarithmic scale. The plots of this intensity profile created with Dubra-Ferrari and BDW are not visually distinguishable; within the telescope aperture, the maximum difference in intensity is 4.8 × 10⁻¹⁵.

Fig. 4 A comparison of the time to run the two propagation algorithms, as a function of number of points along the edge of the entrance aperture.

Download Full Size | PDF

Fig. 5 Intensity at the plane of the telescope aperture following the occulter, plotted on a base-10 logarithmic scale. The telescope aperture is shown by a dashed circle.

Download Full Size | PDF

5. Conclusion

We present a new approach to modeling fields following planet-finding occulters, based on the boundary-diffraction-wave formulation laid out by Miyamoto and Wolf. It has proved to be simple to modify the shape and recalculate with this approach to include various types of occulter errors, as well as fast to evaluate, improving on the current best method by 2–3× with no loss in precision. An efficient implementation has been laid out, and we hope it proves useful to others in the field.

Acknowledgments

This work was carried out as a part of the Imaging Performance Study funded by NASA’s Exoplanet Exploration Program; the research was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

References and links

1. N. J. Kasdin, D. N. Spergel, R. J. Vanderbei, D. Lisman, S. Shaklan, M. Thomson, P. Walkemeyer, V. Bach, E. Oakes, E. Cady, S. Martin, L. Marchen, B. Macintosh, R. E. Rudd, J. Mikula, and D. Lynch, “Advancing technology for starlight suppression via an external occulter,” Proc. SPIE 8151, 81510J (2011). [CrossRef]

2. S. B. Shaklan, M. C. Noecker, T. Glassman, A. S. Lo, P. J. Dumont, N. J. Kasdin, E. J. Cady, R. Vanderbei, and P. R. Lawson, “Error budgeting and tolerancing of starshades for exoplanet detection,” Proc. SPIE 7731, 77312G (2010). [CrossRef]

3. T. Glassman, A. Johnson, A. Lo, D. Dailey, H. Shelton, and J. Vogrin, “Error analysis on the NWO starshade,” Proc. SPIE 7731, 773150 (2010). [CrossRef]

4. R. J. Vanderbei, E. J. Cady, and N. J. Kasdin, “Optimal occulter design for finding extrasolar planets,” Astrophys. J. 665, 794–798 (2007). [CrossRef]

5. E. Cady, L. Pueyo, R. Soummer, and N. J. Kasdin, “Performance of hybrid occulters using apodized pupil Lyot coronagraphy,” Proc. SPIE 7010, 70101X (2008). [CrossRef]

6. C. J. Copi and G. D. Starkman, “The Big Occulting Steerable Satellite [BOSS],” Astrophys. J. 532, 581–592 (2000). [CrossRef]

7. W. Cash, “Detection of earth-like planets around nearby stars using a petal-shaped occulter,” Nature 442, 51–53 (2006). [CrossRef] [PubMed]

8. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996).

9. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999).

10. R. Soummer, L. Pueyo, A. Sivaramakrishnan, and R. J. Vanderbei, “Fast computation of Lyot-style coronagraph propagation,” Opt. Express 15(24), 15935–15951, (2007). [CrossRef] [PubMed]

11. E. Cady, “Design, tolerancing, and experimental verification of occulters for finding extrasolar planets,” PhD thesis, Princeton University, 2010.

12. A. Dubra and J. A. Ferrari, “Diffracted field by an arbitrary aperture,” Am. J. of Phys. 67, 87–92 (1999). [CrossRef]

13. P. Dumont, S. Shaklan, E. Cady, J. Kasdin, and R. Vanderbei, “Analysis of external occulters in the presence of defects,” Proc. SPIE 7440, 744008 (2009). [CrossRef]

14. K. Miyamoto and E. Wolf, “Generalization of the Maggi-Rubinowicz theory of the boundary diffraction wave part I,” J. Opt. Soc. Am 52, 615–622 (1962). [CrossRef]

15. K. Miyamoto and E. Wolf, “Generalization of the Maggi-Rubinowicz theory of the boundary diffraction wave part II,” J. Opt. Soc. Am 52, 626–636 (1962). [CrossRef]

16. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes. The Art of Scientific Computing (Cambridge University Press, 2007).

Boundary diffraction wave integrals for diffraction modeling of external occulters

Abstract

1. Introduction

2. Current occulter modeling techniques

3. The boundary diffraction wave formulation

3.1. Vector potentials for occulter propagation

4. Efficient implementation

4.1. Computational results

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Tables (1)

Equations (49)

Optics Express