## 1.

## Introduction

The Wide-Field Infrared Survey Telescope-Astrophysics Focused Telescope Assets (WFIRST-AFTA) coronagraph is designed to have a powerful high contrast performance for direct imaging of exoplanets by their reflected light. Any coronagraph system considered for the WFIRST-AFTA project should solve the main design problem of removing the starlight diffracted from the telescope central obstruction. The phase-induced amplitude apodization (PIAA) complex mask coronagraph (CMC)^{1} based on lossless pupil apodization by the beam shaping concept^{2} combined with a complex amplitude focal plane mask gives a solution to this problem.

Such a combination can be designed to concentrate the diffracted starlight in the pupil area that matches the telescope obstruction and block it with the Lyot stop that replicates the obstruction shape. It also allows the inner working angle as small as $0.8\lambda /D$ and much milder aspheric mirror shapes in comparison with the original PIAA design. However, it requires a focal plane mask that affects both phase and amplitude, as opposed to the simpler opaque focal plane mask used for PIAA.

The phase-induced amplitude apodization complex mask coronagraph (PIAACMC) (both the mirror shapes and the phase mask) can be optimized in a way that takes into account both chromatic diffraction effects caused by the telescope obstructed aperture, and the tip-tilt sensitivity of the coronagraph,^{3}^{,}^{4} and still attain the contrast of ${10}^{-9}$ between 2.0 and 4.0 $\lambda /D$ with 10% broadband visible light (550 nm)^{4} that satisfies the WFIRST-AFTA requirements.^{5}

The remapping optics of the PIAACMC system can be realized with lenses or mirrors. While either option is suitable at the ${10}^{-7}$ raw contrast level of ground-based systems, mirrors are the preferred option for reaching ${10}^{-8}$ contrast and beyond, as they are free from ghosts and wavefront errors induced by substrate inhomogeneities. Lens-based PIAACMC systems offer the convenience of being on-axis and can be inserted in an existing beam.^{6} Mirror-based systems can only be built on-axis for beams with large central obstruction, in a Cassegrain-like optical configuration. For all other configurations (unobstructed pupil, small central obstruction), mirror-based PIAACMC have to be built off-axis.

The typical design process for PIAA mirrors is to first optimize their shapes for an abstract on-axis collimated space configuration and then convert this design to a real unobstructed off-axis physical configuration like WFIRST-AFTA PIAACMC configuration (see Fig. 1, Ref. 3). The pupil asymmetry introduced by the telescope obstruction produces the system asymmetry that could be shared between the PIAA mirror shapes and the occulter (see the difference between Gen 2 and 3 designs in Ref. 3). Due to this asymmetry and Fresnel propagation effects, the procedure used to design the previous generation of off-axis PIAA optics is not suitable to design the off-axis system considered in this paper.

The goal of this paper is to present a method to calculate PIAA optics for an arbitrary off-axis geometry of the optical system based on the given on-axis optics shapes. By “arbitrary,” we mean different off-axis geometries for the system entrance and exit beams as well as a possible beam compression that the system should perform (one of the requirements for the proposed WFIRST-AFTA PIAACMC configuration^{3}). Though diffraction critically affects the final system performance, our consideration will be based on geometrical optics only, assuming that the diffraction propagation effects are taken into account during the on-axis system optimization and are unaffected by conversion to off-axis. On the contrary to the full diffraction-based approach, the considered procedure is much less calculation intensive and provides the global convergence to the solution. The obtained solution can be used as a starting point for the full diffraction-based approach.

The phase mask design/characterization is outside of the scope of this paper and is considered in a paper by Kern et al.^{3}

## 2.

## Optical Layout

The proposed optical layout of the PIAA system for the WFIRST-AFTA PIAACMC (Fig. 1) includes two PIAA mirrors M1 and M2. In this collimated-beam-to-focus system, the collimated input beam is reshaped by the M1 mirror and is focused by the M2 mirror. The angle of the input beam $\alpha $ is equal to 20 deg (approx.) and the output beam angle $\beta $ is equal to 25 deg (approx.). The input beam diameter is about 40 mm, the distance $d$ between mirrors is about 750 mm, and the focal length $f$ of M2 (the distance between the M2 center and the system focus) is about 1200 mm. The M1 mirror also performs the compression of the input beam needed to match the focal plane point spread function (PSF) size with the diameter of the focal plane occulter. As a result the beam diameter on M2 is approximately equal to 15 mm.

We define a coordinate reference frame shown in Fig. 1 with the origin at the center of PIAA M1 and the $z$-axis parallel to the central ray between M1 and M2 mirrors. All mirror shapes will be expressed in this reference frame as functions of $(x,y)$.

The base shape of M1 and M2 mirrors (i.e., without any PIAA remapping) is a plane and an off-axis parabola (OAP), respectively. All the following mirror shapes are calculated as corrections to these shapes that are given by

## (1)

$$\mathrm{PLANE}(x,y)=y\times \mathrm{tan}(\alpha /2),\phantom{\rule{0ex}{0ex}}\mathrm{OAP}(x,y)=\frac{{f}^{2}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\beta -{x}^{2}-{(y-f\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\beta )}^{2}}{2f(1+\mathrm{cos}\text{\hspace{0.17em}}\beta )}+d.$$Two additional mirrors are used to perform remapping calculations. They are the flat mirror FLAT and the off-axis parabola ${\mathrm{OAP}}_{f}$

## (2)

$$\mathrm{FLAT}(x,y)=\frac{2d}{\mathrm{cos}\text{\hspace{0.17em}}\alpha}+(y+\frac{2d}{\mathrm{sin}\text{\hspace{0.17em}}\alpha})\mathrm{tan}(\alpha /2),\phantom{\rule{0ex}{0ex}}{\mathrm{OAP}}_{f}(x,y)=\frac{-{f}^{2}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\beta +{x}^{2}+{(y-f\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\beta )}^{2}}{2f(1+\mathrm{cos}\text{\hspace{0.17em}}\beta )}+d-2f\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\beta .$$## 3.

## Off-Axis Mirror Shapes

In this section the procedure to calculate off-axis PIAA shapes is given. The shape of both mirrors can be presented as a sum of three terms

## (3)

$${M}_{1}(x,y)=\mathrm{PLANE}(x,y)+{f}_{1}(r)+{\mathrm{OAT}}_{1}(x,y),\phantom{\rule{0ex}{0ex}}{M}_{2}(x,y)=\mathrm{OAP}(x,y)+{f}_{2}(r)+{\mathrm{OAT}}_{2}(x,y).$$## 3.1.

### Remapping Functions

The ray-tracing calculations of the PIAA optics are based on the remapping function concept. For a source on the optical axis the remapping function establishes a correspondence between a ray position ${\mathbf{r}}_{\mathbf{1}}=({x}_{1},{y}_{1})$ in the entrance pupil and the position ${\mathbf{r}}_{\mathbf{2}}=({x}_{2},{y}_{2})$ of the same ray in the exit pupil

Note, that the origin of the reference frames in Eq. (4) matches with the position of the central ray in the entrance and exit pupils.Since we take the on-axis mirror shapes as given, the ray position in the entrance and exit pupils, and thus two-dimensional (2-D) remapping functions ${x}_{2}({x}_{1},{y}_{1})$ and ${y}_{2}({x}_{1},{y}_{1})$, can be readily calculated with a ray-tracing procedure.

The starting point to design mirror shapes for the WFIRST-AFTA coronagraph was the diffraction-based optimized solution (Fig. 2, Ref. 3; “revised” PIAACMC design, Ref. 4) obtained for an on-axis collimated-beam-to-collimated-beam PIAA system without beam compression. Although this solution is not circularly symmetric because of the central obstruction asymmetry, the deviation from the symmetry for both mirrors is small. That makes it possible to define the azimuthally averaged remapping function

## (5)

$${r}_{2}({r}_{1})={\langle \sqrt{{x}_{2}^{2}({r}_{1}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\varphi ,{r}_{1}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\varphi )+{y}_{2}^{2}({r}_{1}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\varphi ,{r}_{1}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\varphi )}\rangle}_{\varphi},$$The remapping function ${r}_{2}({r}_{1})$ constructed in according with Eq. (5) provides a good initial approximation to calculate the off-axis mirror shapes for all the rays within the working beam. However, the discontinuity of derivatives of the remapping function at the beam edge and the absence of the remapped values outside of the working beam creates numerous difficulties for the ray-tracing calculations that are suppose to work with real square arrays of points.

To avoid these problems we created an extended remapping function ${R}_{2}({r}_{1})$ (Fig. 2) by transposing the remapping function ${r}_{2}({r}_{1})$ onto the intervals $[{r}_{\mathrm{b}1},R]$ and $[{r}_{\mathrm{b}2},R]$ backward, beyond the working beam radii ${r}_{\mathrm{b}1}$ (on PIAA M1) and ${r}_{\mathrm{b}2}$ (on PIAA M2), and making sure that the derivatives match on the beam edge. The ${R}_{2}({r}_{1})$ is given as follows:

## (6)

$${R}_{2}({r}_{1})={r}_{2}({r}_{1})\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{for}\text{\hspace{0.17em}\hspace{0.17em}}{r}_{1}<{r}_{\mathrm{b}1},\phantom{\rule{0ex}{0ex}}\{\begin{array}{c}{r}_{1}={r}_{\mathrm{b}1}+(1+{\alpha}_{1}{r}_{1}^{\prime}){r}_{1}^{\prime}\\ {r}_{2}={r}_{\mathrm{b}2}+(1+{\alpha}_{2}{r}_{2}^{\prime}){r}_{2}^{\prime}\\ {r}_{\mathrm{b}2}-{r}_{2}^{\prime}={R}_{2}({r}_{\mathrm{b}1}-{r}_{1}^{\prime})\end{array}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{for}\text{\hspace{0.17em}\hspace{0.17em}}{r}_{\mathrm{b}1}\le {r}_{1}\le R,$$It should be noted that for a system with beam compression ${r}_{\mathrm{b}2}$ is not equal to ${r}_{\mathrm{b}1}$. The desired beam compression can also be included in the final expressions for 2-D remapping functions ${X}_{2}({x}_{1},{y}_{1})$ and ${Y}_{2}({x}_{1},{y}_{1})$ through the rescaling factor ${r}_{\mathrm{b}2}/{r}_{\mathrm{b}1}$

## (7)

$${X}_{2}({x}_{1},{y}_{1})=\frac{{r}_{\mathrm{b}2}}{{r}_{\mathrm{b}1}}{x}_{2}({x}_{1},{y}_{1}),\phantom{\rule{0ex}{0ex}}{Y}_{2}({x}_{1},{y}_{1})=\frac{{r}_{\mathrm{b}2}}{{r}_{\mathrm{b}1}}{y}_{2}({x}_{1},{y}_{1}),$$## 3.2.

### Optics Shapes Calculation

The following iteration procedure is proposed to calculate off-axis PIAA mirror shapes given on-axis ones:

1. Calculate the 2-D remapping functions ${x}_{2}({x}_{1},{y}_{1})$ and ${y}_{2}({x}_{1},{y}_{1})$ by ray-tracing on-axis mirror shapes and azimuthally average them to get the radial remapping function ${R}_{2}({r}_{1})$.

2. Calculate the 2-D remapping functions ${X}_{2}({x}_{1},{y}_{1})$ and ${Y}_{2}({x}_{1},{y}_{1})$ for the compressed beam.

3. Determine base shapes $\mathrm{PLANE}(x,y)$ and $\mathrm{OAP}(x,y)$ of mirrors per Eq. (1).

4. Calculate radial mirror profiles ${f}_{1}(r)$ and ${f}_{2}(r)$ for the on-axis collimated-beam-to-focus system with the radial remapping ${R}_{2}({r}_{1})$ (Sec. 3.3).

5. Start a loop. For $i$’th iteration step construct an estimate ${\mathrm{OAT}}_{1,i}(x,y)$ as a sum of small corrective terms $d{M}_{1,j}(x,y)$ calculated during previous iteration steps (Sec. 3.4)

6. Determine the off-axis PIAA term ${\mathrm{OAT}}_{2,i}(x,y)$ that provides equal optical path length (OPL) for all the rays propagated from the source to the main PIAA focus. This step is performed assuming that the off-axis PIAA term ${\mathrm{OAT}}_{1}(x,y)$ is equal to ${\mathrm{OAT}}_{1,i}(x,y)$.

7. Close the loop by returning to step 5.

## 3.3.

### Radial Mirror Terms

For the case of a radially symmetric collimated-beam-to-collimated-beam PIAA system, the mirror shapes can be directly computed from ${R}_{2}({r}_{1})$ by numerically integrating a pair of differential equations.^{2}

In this section, we consider the case of the on-axis collimated-beam-to-focus radially symmetric system. The proposed solution is similar to the procedure we use to calculate the 2-D off-axis correction to the M1 shape (Sec. 3.4).

The radial mirror profiles ${f}_{1}(r)$ can be derived from

if the derivative $\mathrm{d}{f}_{1}(r)/\mathrm{d}r$ is known. The derivative $\mathrm{d}{f}_{1}(r)/\mathrm{d}r$ can be easily estimated for a collimated-beam-to-focus system by using the “desired” normal [i.e., one that generates the desired remapping in according with Eq. (5)] to the M1 surface.In according with the reflection law, the “desired” normal ${\mathbf{n}}_{\mathbf{i}}=({n}_{r}^{i},{n}_{z}^{i})$ to the surface of the M1 mirror can be written as

## (10)

$${\mathbf{n}}_{\mathbf{i}}{(\mathbf{r}}_{\mathbf{1}})=({n}_{r}^{i},{n}_{z}^{i})=\frac{{\mathbf{r}}_{\mathbf{1}\mathbf{s}}}{{r}_{1s}}+\frac{{\mathbf{r}}_{\mathbf{12}}}{{r}_{12}},$$## 3.4.

### Calculation of $\mathrm{d}{M}_{1}(x,y)$

To calculate the corrective term $\mathrm{d}{M}_{1,j}(x,y)$ (Sec. 3.2) we should

1. For given 2-D remapping functions [Eq. (7)] determine the “desired” normal ${\mathbf{n}}_{\mathbf{i}}=({n}_{x}^{i},{n}_{y}^{i},{n}_{z}^{i})$ to the M1 surface assuming that M1 and M2 mirror shapes match the current shapes approximation. Similarly Eq. (10), the “desired” normal

where vectors ${\mathbf{r}}_{\mathbf{1}\mathbf{s}}$ and ${\mathbf{r}}_{\mathbf{12}}$ are determined in Sec. 3.3. To determine a particular ray position on the M2 surface, the geometrical propagation in direction from the exit pupil to the off-axis parabola ${\mathrm{OAP}}_{f}$ and then to the M2 mirror should be used. The ray position in the exit pupil can be determined as the remapped [in accordance with Eq. (7)] ray position in the entrance pupil. By the entrance pupil, we mean a plane upstream of the flat mirror $\mathrm{FLAT}(x,y)$ [Eq. (2)]. The exit pupil is the plane downstream of the off-axis parabola ${\mathrm{OAP}}_{f}(x,y)$ [Eq. (2)] optically conjugated with the M2 mirror. Both the entrance and exit pupils are orthogonal to the collimated input/output beams and all the rays are considered to be collimated between the entrance pupil and the flat mirror, and between the exit pupil and ${\mathrm{OAP}}_{f}$.## (12)

$${\mathbf{n}}_{\mathbf{i}}(\mathbf{x},\mathbf{y})=({n}_{x}^{i},{n}_{y}^{i},{n}_{z}^{i})=\frac{{\mathbf{r}}_{\mathbf{1}\mathbf{s}}}{{r}_{1s}}+\frac{{\mathbf{r}}_{\mathbf{12}}}{{r}_{12}},$$2. Obtain the “current” normal (i.e., gradient vector of the surface) ${\mathbf{n}}_{\mathbf{c}}=({n}_{x}^{c},{n}_{y}^{c},{n}_{z}^{c})$ to the current M1 mirror shape.

3. Calculate derivatives $\partial \mathrm{d}{M}_{1,i}(x,y)/\partial x$ and $\partial \mathrm{d}{M}_{1,i}(x,y)/\partial y$

4. Calculate $\mathrm{d}{M}_{1,i}(x,y)$ through a line integral

## (15)

$$\mathrm{d}{M}_{1,i}(x,y)={\int}_{(\mathrm{0,0})}^{(x,y)}\frac{\partial \mathrm{d}{M}_{1,i}(x,y)}{\partial x}\mathrm{d}x+\frac{\partial \mathrm{d}{M}_{1,i}(x,y)}{\partial y}\mathrm{d}y.$$Since, for any physically possible remapping, the integration in Eq. (15) along any path that connects the mirror center with the point $(x,y)$ should give the same result, $\mathrm{d}{M}_{1,i}(x,y)$ values can also be averaged for different integration paths.

5. To reduce the M1 errors caused by the numerical integration/interpolation, an appropriate smoothing algorithm should be applied for each $\mathrm{d}{M}_{1,i}(x,y)$ estimation.

## 3.5.

### Calculation of OAT_{2} (x, y)

To determine ${\mathrm{OAT}}_{2,i}(x,y)$, we use an iterative procedure II that corrects the OPL of the output rays for every step of the main iteration. Note that iteration procedure II is different from the main iteration discussed earlier in Sec. 3.2.

During $k$'th step of the iteration II, the function $\mathrm{d}{M}_{2,k}({x}_{2},{y}_{2})$ is calculated

## (16)

$$\mathrm{d}{M}_{2,k}({x}_{2},{y}_{2})=-\frac{{\mathrm{OPL}}_{k}({x}_{2},{y}_{2})-\mathrm{OPL}(0,0)}{2},$$## (17)

$${\mathrm{OAT}}_{2,i}({x}_{2},{y}_{2})=\sum _{k}\mathrm{d}{M}_{2,k}({x}_{2},{y}_{2}),\phantom{\rule{0ex}{0ex}}{\mathrm{OAT}}_{2,0}({x}_{2},{y}_{2})=0.$$One question, related to the iterations described above needs to be mentioned specially. During each iteration step numerous derivative calculations are performed to estimate the normal to the mirror surfaces at the location of each ray. The integration, differentiation, and interpolation executed during each iteration step produce and propagate numerical errors. Those errors result in so-called numerical catastrophe when the calculated derivatives show infinite growth starting from some iteration number. To avoid the numerical catastrophe, the data arrays related to calculation of the $\mathrm{d}{M}_{1}(x,y)$ and $\mathrm{d}{M}_{2}(x,y)$ terms should be smoothed with any appropriate algorithm. It is clear though that the convergence of iterations should slow down (or even be broken), if the spatial scale of the smoothing becomes comparable with the expected size of features on the mirror surfaces. For the system considered here, the minimal mirror features size is determined by the diffraction-like structure described in Sec. 4. As a result, the residual output beam OPD is constrained to the ${10}^{-10}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ level, which is better than the physical limit of optics manufacturing. Also, note that the ${10}^{-10}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ RMS errors in mirror shapes produce 40 to $50\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ ray scattering in the focal plane of the system, i.e., comparable with the size of the diffractive PSF of the system (50 to $60\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ depending on the wavelength) and a $50\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ beam walking in the M2 mirror plane.

## 4.

## Mirror Shapes

The calculated mirror shapes for the WFIRST-AFTA coronagraph are shown in Fig. 4. To present the deviation of shapes from common optical shapes we presubtracted following terms

## (18)

$${P}_{1}(x,y)=0.171215209\times y+0.24227703\times ({x}^{2}+{y}^{2})\phantom{\rule[-0.0ex]{2em}{0.0ex}}\mathrm{a}\mathrm{n}\mathrm{d}\phantom{\rule{0ex}{0ex}}{P}_{2}(x,y)=0.228083501\times y+0.64898643\times ({x}^{2}+{y}^{2})$$An additional point to emphasize is that the diffraction propagation introduces radial oscillations both in amplitude and phase of the output beam (Figs. 5 and 6) which the initially optimized on-axis mirror shapes take into account and correct for.

As a result, the geometrical OPD of the output beam is not equal to zero as expected in according to the pure geometrical optics model (Fig. 5), although it appears to be equal to zero if we account for diffraction propagation. To take into account the above effect, an additional term $\mathrm{PHASE}(x,y)$ for the M2 mirror shape should be introduced

## (19)

$${M}_{2}(x,y)=\mathrm{OAP}(x,y)+{f}_{2}(r)+{\mathrm{OAT}}_{2}(x,y)+\mathrm{PHASE}(x,y)\mathrm{.}$$## (20)

$$\mathrm{d}{M}_{2,k}({x}_{2},{y}_{2})=-\frac{{\mathrm{OPL}}_{k}({x}_{2},{y}_{2})-\mathrm{OPL}(0,0)-{\mathrm{OPD}}_{d}(x,y)}{2},$$## 5.

## Ray-Tracing Check

The derived optical shapes with the “phase” term $\mathrm{PHASE}(x,y)=0$ have been checked with the ray-tracing code that has been developed to simulate the previous generations of PIAA optics manufactured by AXSYS Inc.^{7} and Tinsley Inc.^{8} Results of the ray-tracing simulation are shown in Fig. 6 and include the wavefront amplitude and phase in the exit pupil of the system calculated in the geometrical optics approximation.

Though the wavefront amplitude looks noisy with a 4.7% RMS white noise in amplitude caused by the limited number of rays used for the simulation (about $115\text{\hspace{0.17em}\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rays}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{per}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{pixel}}^{2}$), the calculated values show close correspondence with the designed exit pupil amplitudes except for the area that is blocked by the telescope central obstruction. At the same time, the expected profile seems to be deformed by the clearly observable diffraction-like structure. The presence of this structure [as well as the presence of the “phase” term $\mathrm{PHASE}(x,y)$] is easy to understand if one takes into account the diffraction-based optimization of the on-axis mirror shapes mentioned above. This optimization produces the shapes that deliver the designed “smooth” (without diffraction-like structure) amplitude and phase profiles in the presence of the Fresnel diffraction. As a result, the “smooth” shape of the mirrors as well as the exit pupil amplitude and phase profiles (in the geometrical optics approximation) appears to be modulated with the “inverse” (in respect to the Fresnel diffraction) diffraction-like structure that is responsible for the correction of diffraction propagation effects.

In an absence of the $\mathrm{PHASE}(x,y)$ term, the exit pupil phase appears to be flat with the wavefront errors of about $\lambda /3000$. The exit pupil phase errors are only limited by the residual output beam OPD with a ${10}^{-10}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ RMS discussed in Sec. 3.5. The value of ${10}^{-10}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ seems to be a reasonable estimate for the accuracy of the PIAA mirror shapes calculation with the described algorithm.

## 6.

## Summary

In this paper, we presented a method to calculate off-axis mirror shapes for the PIAACMC assuming that the on-axis mirror shapes are known. The method is based on geometrical ray tracing and is able to calculate off-axis PIAA mirror shapes for an arbitrary geometry of the input and output beams. The accuracy of off-axis PIAA shapes calculation is limited to the value of ${10}^{-10}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, which is better than the physical limit of optics manufacturing.

## Acknowledgments

This work was supported by the National Aeronautics and Space Administration’s WFIRST-AFTA program. It was carried out at the NASA Ames Research Center in collaboration with NASA Jet Propulsion Laboratory. The authors would like to thank Brian Kern for precursor work and helpful discussions. The preliminary version of this paper has been published in Proceedings of SPIE 9605, 2015.^{9}