Implementation and benchmarking of a crosstalk-free method for wavefront Zernike coefficients reconstruction using Shack-Hartmann sensor data

In wavefront characterization, often the combination of a Shack-Hartmann sensor and a reconstruction method utilizing the Cartesian derivatives of Zernike circle polynomials (the least-squares method, to be called here Method A) is used, which is known to introduce crosstalk. In \citep{janssen2014zernike} a crosstalk-free analytic expression of the LMS estimator of the wavefront Zernike coefficients is given in terms of wavefront partial derivatives (leading to what we call Method B). Here, we show an implementation of this analytic result where the derivative data are obtained using the Shack-Hartmann sensor and compare it with the conventional least-squares method.


Introduction
In applications involving wavefront reconstruction, Shack-Hartmann sensors are frequently used. An often-used polynomial basis to expand the phase of the scalar field in an optical system is the set of Zernike circle polynomials. The wavefront Zernike coefficients can be estimated with for example, the least-squares method, to be called Method A henceforth. However, already since the 1980s, it is known that there is a fundamental problem in the combination of the least-squares method and Zernike polynomials: this method presents crosstalk between coefficients. In 2014, Janssen introduced an analytic result in the form of a new relation between the wavefront derivative data and the wavefront Zernike coefficients (Janssen, 2014). This relation theoretically solves the crosstalk problem that is present in the least-squares method. In particular, the wavefront coefficients, once estimated, do not change anymore when the number of Zernike circle polynomials involved in the fit is increased further.
In this paper, we show an implementation of Janssen's analytic result in the reconstruction of several wavefronts that have been generated using a spatial light modulator. In particular, we show that this method, to be called Method B henceforth, does not suffer from cross-talk between the coefficients.
The paper is organized as follows. In Section 2 we present theoretical aspects of the problem, and in Section 3, the experimental setup and methods are explained. The actual results are presented and discussed in Section 4, and the conclusions from these measurements are presented in Section 5.

Theory
In this section, Zernike circle polynomials are discussed and used as a basis to describe the wavefront deviation with respect to a suitably defined reference sphere of an optical system. This wavefront deviation or wavefront aberration function will be loosely called 'wavefront aberration' or 'aberration' in what follows. We start this section by showing the definition of Zernike polynomials and a description of the data that is obtained with the Shack-Hartmann sensor is discussed. Next, we discuss two reconstruction methods: the well-known least-squares method (Method A) and the recently method introduced by Janssen(Janssen, 2014) (Method B).

Zernike polynomials as a basis for wavefront expansion
In systems where aberrations are desired to be known, such as microscopes and telescopes, a circular aperture is usually present. In order to describe these aberrations, Zernike circle polynomials are commonly used. In this paper, the American National Standards Institute (ANSI) definition for Zernike circle polynomials is considered, which is commonly used to describe wavefront reconstruction from Shack-Hartmann data. 1 As defined in Thibos et al. (2002), the polynomials in polar coordinates (ρ, θ) are given by such that z m n (ρ, θ) is a real-valued, orthonormal (on the unit disc) expression for Zernike circle polynomials, n = 0, 1, 2, ... is the degree of the Zernike polynomial, and m = 0, ±1, ±2, ... the azimuthal order satisfying: (2.5) n − |m| is even, (2.6) |m| ≤ n, (2.7) ρ ≤ 1 being the radius on the unit disc, and δ nn the Kronecker delta function. Here, the n th "order" polynomial is referred to as the n th degree polynomial.
As mentioned above, these circle polynomials are orthonormal on the unit disc, i.e., 1 π 1 0 2π 0 z m n (ρ, θ)z m n (ρ, θ)ρdθdρ = δ nn δ mm . (2.8) With this definition, any real-valued wavefront function (defined on the unit disc) can be described as a linear combination of the Zernike circle polynomials given by where η m is the set of allowed n values dependent on m, namely η m = {|m|, |m| + 2, |m| + 4, . . .}, when m = 0, and η m = {2, 4, . . . , ∞} when m = 0 and z 0 0 = 1. This set ensures that the constraints set in n and m are all met. Finally, a m n is the real-valued Zernike coefficient.

Complex-valued definition of the Zernike circle polynomials
Following the complex-valued definition of the Zernike circle polynomials as used in (Janssen, 2014) in polar coordinates (ρ, θ) on the unit disc, we have that (2.10) The radial polynomial R |m| n is given as in Equation 2.1 above, and can also be given as is the Jacobi polynomial of degree k, which is orthogonal with respect to the weight (1 − x) α (1 + x) β on the interval [−1, 1]. Note that the factor ρ |m| is missing in (Janssen, 2014), Eq. (6), but it is merely a typo and has no consequences in the further developments in that reference. We set Z m n = 0 for all values of n and m where n − |m| is odd or negative. There is the normalization condition This orthogonality means that, like the real-valued Zernike polynomial, any sufficiently smooth (complex) wavefront can be described by where α m n are generally complex-valued coefficients corresponding to the complex Zernike polynomial and η n is as in Equation 2.9. Expanding the complex exponential in Equation 2.10 leads to the following conversion between this definition of complex Zernike polynomials and the ANSI standard (2.14) where N m n is defined in Equation 2.2. From the definition of the complex Zernike polynomial in Equation 2.10 we can also see that the complex conjugate Z |m| n * is equal to Z −|m| n . This observation, together with the relations between the real and complex Zernike polynomials leads to the expression (2.15) When using complex polynomials to describe a real wavefront, the coefficients are usually complex valued. In order to transform them back to the coefficients of the real-valued Zernike z m n , the following relations can be used (2.16) Here, the Zernike polynomials are evaluated on the computer, and data is saved in vectors and matrices. To make it easy to loop over all polynomials, a single index is introduced. In the notation by Thibos et al. (2002), this single index j = 0, 1, ... is given by where n is the Zernike degree and m the azimuthal order. From now on, we will write Z j = z m n . The reverse can also be done, that is finding n and m from j as follows where x denotes the "ceiling" function, that is the smallest integer greater or equal than x.

Data from the Shack-Hartmann sensor
It is well known that aberrations in an optical system degrade its imaging quality. The phase aberration is defined by a function ϕ W = arg(W ). A typical way of measuring this function is by using a Shack-Hartmann sensor. A Shack-Hartmann sensor consists of a camera chip and lenslet array. The lenslets are placed at the focal distance from the camera chip such that an incoming plane wave will be focussed as many spots on the camera. If the incoming wave is aberrated, the position of each spot on the sensor will changed. It can be derived from the theory how much a spot is displaced due to an aberration. Because this displacement is linear, the reverse problem can be solved as well. In other words, one can retrieve the aberration given the spot displacement. Dai (2008) has proven that for a varying wavefront over the subaperture, the slope needs to be averaged over that sub-aperture. The expressions for the average slopes in terms of displacement of the spot on the camera becomes where ∆x, ∆y are the shift in the x− and y− positions of the spots, r is the radius of the incoming beam on the Shack-Hartmann sensor, and f is the focal length of the sensor lens array, Σ is the illuminated sub-aperture domain, with surface area A Σ . Here Σ and A Σ change when the sub-aperture is only partially illuminated (i.e., at the edge of the beam). This averaging of the slope is important in recovering the wavefront.

Method A
The least-squares (LSQ) fit is based on the real Zernike polynomials, and uses the fact that the coefficients are not dependent on x and y. The wavefront in x, y coordinates (with x 2 + y 2 ≤ 1) can be described as a m n z m n (x, y). (2.21) The measured quantity, however, is not W but ∂W ∂x and ∂W ∂y . Taking the partial derivatives to x and y results in the over determined system of (2.22) This system can be solved for a finite set of polynomials. Using Equation 2.20, given the displacements ∆x and ∆y, one can create a vector s containing the slopes as such (2.23) The partial derivatives of the Zernikes in the x-and y-direction can also be put in a matrix, called the geometry matrix. We recall the convention that Z j = z m n with j, m, n related as in Equation 2. .20. It should be noted that the positions over which the averaging is done is normalized to the unit disc. These windows are the same defined in Equation 2.20. The matrix will have a size of (2n spot × J), where J is the maximum index for the Zernike modes used to have a good approximation of the true wavefront. The expression for G becomes (2.24) The system of Equation 2.22 can then be written as where a is the vector containing the Zernike coefficients. The least-squares estimation of a m n becomes a ≈ G + · s, (2.26) where G + the generalized inverse of the geometry matrix. This is an approximation as G only contains the information of a finite number of Zernike modes, and their contribution is averaged over the lenslet array.

Method B
Method B relies on an analytical relation found between the local derivatives of the wavefront and Zernike polynomials. This is in contrast with method A, where there is a link between the local derivatives of the wavefront and the derivatives of the Zernike polynomials. The reconstruction using Method B is based on the identities: where these identities allow the expression of any derivative Zernike polynomial as a sum of Zernike polynomials. Here we use the convention of Equation 2.10 and Equation 2.11 for the Zernike circle polynomials. (Janssen, 2014) has found that the LMS complex coefficients are given as (see appendix):α and where δ nn is the Kronecker delta equal to 1 if n = n and 0 otherwise. Furthermore β + and β − are the Zernike coefficients of ∂W ∂x ± i ∂W ∂y so that Note that the LMS coefficientsα m n are analytically related to four β coefficients, namely (β + ) m+1 n+1 , (β + ) m+1 n−1 , (β − ) m−1 n−1 , and (β − ) m−1 n+1 . As a consequence of orthogonality of the Z m n in the right-hand side of Equation 2.31, the LMS coefficients, once estimated, do not change anymore when the number of Zernike terms used in a finitized version of Equation 2.31 is increased further.
The fact thatα m n is analytically related to β coefficients is desirable, because β-coefficients can directly be estimated (in a least-squares sense) from measurable quantities. This fit is made in the same way as the leastsquares method. This means that also the complex Zernike polynomials need to be averaged over the lenslets. The fit to get the β coefficients, however, is done with a different basis than in the least-squares method to find the a-coefficients. The effects of this are discussed in the following section.
As a note for this method, when n = |m|, there will be non-existent combinations of n and m in Equation 2.30. In that case the value of β will be set to 0. For instance,α 1 1 is among others dependent on (β + ) 2 0 , which goes against the constraint |m| ≤ n. To go from complex coefficient α m n to the real coefficient a m n , the relations in Equation 2.16 can be used. In (Mahajan and Acosta, 2017), a seemingly different approach, based on vector polynomials, is used to express the wavefront Zernike coefficients in terms of wavefront slope data. However, we show in the appendix that this method is essentially equivalent to Janssen's algorithm in (Janssen, 2014).

The main difference between Method A and Method B
The main difference between Method A and Method B for finding the coefficients is in how the fitting is done. Both methods use a least-squares fit using a geometry matrix, but the matrix elements are constructed differently. In Method A, the geometry matrix elements are evaluations of the average gradient of the real-valued Zernike polynomials over certain subdisks, while for Method B it is the average of the complex-valued Zernike polynomials over the same windows.
The gradients of the Zernike polynomials are known not to be orthogonal. This can cause problems called crosstalk when fitting the coefficients, especially when there are more aberrations present in the system than are being fit.
If a is an M -dimensional vector containing the coefficients of the aberrations present in the system, the slopes on the Shack-Hartmann sensor can be determined as s = Ga, where s is an 2n spot long vector containing the x-and y-displacement on the Shack-Hartmann sensor and G an 2n spot × J geometry matrix defined in Equation 2.24.
The notation Z j = z m n is used since j, n, m are related (see Eq. Equation 2.17-Equation 2.19). Z 0 is not included in the matrix since Z 0 = z 0 0 = 1, and its partial derivatives are equal to 0. Note that the first column contains the x-and y-derivatives of the first Zernike polynomial evaluated in all n spot points. When a least-squares estimation of the coefficientsâ (where the hat means to indicate that it is an estimated parameter) is made using less Zernike polynomials, up to Zernike polynomial J < M , crosstalk will occur. The estimatorâ can be expressed aŝ where G l is the geometry matrix containing the columns of the first J Zernike polynomials. The estimator will estimate the lower-order values of the coefficients with influence of the higher-order values, because the matrix G + l G will not be an identity matrix.
When estimating the coefficient a m n , a higher-order aberration a m n will influence the estimation if it is not accounted for in G l (i.e., the single index of a m n j > J) and if where in both cases η m and η m are the sets of allowed values for n and n dependent on m and m such that Equations 2.5, 2.6 and 2.7 are all met. Because of this dependence of an expansion coefficient on the maximum degree of the system of equations, it is expected that methods such as the least-squares method will incorrectly estimate the coefficients when there are higher-order aberrations present that are not accounted for in the geometry matrix G l . For Method B, the geometry matrix contains the Zernike polynomials themselves, and therefore it is not expected to present any crosstalk. This is experimentally verified and shown in this paper.

Experimental setup and methods
In Figure 3.1, a schematic view of the experimental setup is shown. The phase of an uniform, collimated laser beam (HeNe laser) is modified by a spatial light modulator (Holoeye PLUTO-2-VIS-056) in such a way that known aberrations are added to the beam. In order to remove the zerothorder reflected light from the SLM, a phase ramp is added to all SLM patterns so that the unaberrated spot is blocked by the iris after being focussed by lens L3. The lenses L1 and L2 and L3 and L4 form two 4f systems.
Using the described experimental setup, all steps necessary to implement method A and Method B can be performed. These are shown in Algorithm 1.
Algorithm 1 Complete measurement and comparison of Shack-Hartmann phase retrieval algorithms 1: Remove initial aberrations of the entire system 2: Add controlled aberration a to the SLM 3: Gather flat and aberrated wavefront Hartmannograms, i.e., the image of the spot pattern generated by the Shack-Hartmann sensor 4: Find the optimal center and radius position for both methods separately such that the RMS value is minimized 5: Calculate Zernike coefficients using optimized center and radius with a desired maximum degree 6: Calculate error between reference and estimate wavefront 3.1 Addressing a phase pattern to the SLM In order to be sure that only the aberration added on the SLM is measured, first all initial aberrations are eliminated. These aberrations can include alignment errors and the surface of the SLM itself, which might not be completely flat (Matsumoto et al., 2008). Removing the initial aberrations has been done by using the Shack-Hartmann sensor and retrieving the wavefront using the least-squares method. The way a phase pattern is addressed to the SLM is as follows. First, an aperture on the SLM is defined. In the current research this is a circular aperture, but the same method is valid if an annular aperture is used. All pixels i within this aperture are used. A "Zernike matrix" can be set up, such that for each pixel i within the aperture the value of all necessary Zernike polynomials can be computed. In matrix form this would be where J is the total number of Zernike polynomials evaluated, and I is the total number of pixels within the aperture. Due to the cyclic nature of the phase pattern and the limits of the SLM, the phase difference assigned to the SLM should be between 0 and 2π. This can be done using the modulo (or mod) operation: a mod n = a − n a/n . If p is a vector containing the values of the individual pixels of the SLM, it can be constructed from the vector a containing the coefficients for the to-be-added aberration by p = (Za) mod 2π.
(3.2) Therefore, in the case of updating the SLM pattern we use where γ is the gamma function of the SLM, calibrated for the used wavelength.
After the correction is done and the wavefront from the SLM is flattened, the phase pattern of the "flat" phase is saved. To this phase pattern, a phase ramp is added to separate the 0 th and 1 st order of the SLM together with the desired aberration, in the same fashion as the correction is added. After this, it is necessary to check if the the maximum phase change over 4 pixels is not being exceeded. Such a check is necessary to avoid aliasing of the SLM-phase. In the current research, this aliasing constraint is simplified to the constraint that the difference between two neighboring pixels should not exceed 0.5π. This is evaluated by letting p i,j be the value of the pixel located at position i, j on the SLM electrode matrix. Then, first two matrices are constructed: (3.4) Afterward, the element-wise minimum is taken between ∆p and 2π − |∆p| for both x and y, in order to account for the modulated phase. If any of the values of this piecewise minimum is above 0.5π, it is said to break the aliasing constraint.

Constructing the matrix G
Recalling from the theory described in the previous section, in order to recover the coefficients describing the aberrated wavefront, s and G have to be constructed. Here, s relies on the two Hartmannograms, one of the flat wavefront and one of the aberrated wavefront from the SLM, and the radius r SH of beam hitting the Shack-Hartmann sensor.
To construct G, it is necessary to average the gradients of the Zernike polynomials or their gradients. The window over which the polynomial has to be averaged can be seen as a scaled version of the lenslet, scaled so that all illuminated lenslets fit the unit disc. In order to compute the Zernike polynomials in these windows, the center position c and the radius r SH of the beam on the Shack-Hartmann sensor have to be known. The window size is estimated by the average distance between the nearest neighbor spots on the Shack-Hartmann sensor. The center position c follows the average position of all the midpoints of the spots on the Shack-Hartmann sensor.
The radius is estimated by calculating the length between the center and the furthest spot from the center.
After this first estimation of the center position and radius is made, a better estimation can be find by optimization. For this optimization, an error has been defined that can be minimized. In this research, a root-meansquare, or RMS type error is chosen. Using the known added aberration and the measured aberration coefficients, an RMS error can be defined. To this purpose, a new Zernike matrix similar to Equation 3.1 is constructed, this time with N points on a grid within the unit disc. The reference phase p ref and the recovered phase p rec can be constructed as where for p ref , the reference vector a ref is used, while for p rec the estimated coefficient vectorâ is used. Using this definition for the reference and recovered phase, the RMS error is determined as where x 2 is the Euclidean vector norm of vector x. This RMS error is then minimized for center position and radius of the beam on the Hartmannogram. A limited memory bound Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) minimization algorithm is applied to find the center position and radius (Fletcher, 1987). The termination conditions for this optimization are: max{|proj(g i )| i = 1, ..., n} ≤ 10 −5 , (3.8) where f k is the value of the RMS error of the k th iteration of the minimization algorithm, and proj(g i ) is i th component of the projected gradient where n projections are made. If any of these statements were true, the optimization was terminated. After looking at the RMS error landscapes, it was found that not every minimum found was a global minimum. If this was the case, the global minimum coordinates were estimated using the RMS landscape graphs, and a brute force optimization was run around those coordinates. This brute force optimization calculates the RMS error value in a grid of points. From the coordinates with the lowest RMS error, a new downhill simplex minimization algorithm is started. This way the global minimum was attempted to be found, and the optimal center and radius positions were determined. In the current research, the necessary aberration coefficients are obtained using a set of Zernike polynomials with a maximum degree of eight. From Soloviev and Vdovin (2006) it is concluded that a good rule of thumb for the amount of Zernike polynomials that can be fit given k spots on the Shack-Hartmann sensor is k /3. In the current research, there were about 144 spots on the Shack-Hartmann sensor, and there are 45 polynomials in the Zernike expansion with a maximum degree of eight.
For method A, the amount of polynomials fit is equal to the amount of retrieved coefficients, while for Method B, it is not so, since for the latter if a fit is made with up to maximum degree of eight, the first seven degrees can be retrieved.
It should be noted that this optimization can take long due to the fact that the geometry matrix needs to be calculated in every iteration, as the values in the matrix depend on the center position and radius.
These optimized parameters can then be used to find the RMS value fitting any degree of Zernike polynomials, and can also be used to determine the error landscape by calculating the RMS values when the center and radius differ slightly from the optimal value.

Results
The two methods have been implemented for three different test cases: a) single Zernike aberrations, b) a combination of three aberrations, assumed not to show crosstalk, and c) four cases of aberrations where crosstalk is present. The aberrations are listed in Table 4.1.     Based on the errors in the reconstruction found in Figure 4.1, Figure 4.2 and Figure 4.3, one can see that both methods provide comparable accuracy in reconstructing the wavefront. However, it is known that Method A presents crosstalk of coefficients when less Zernike powers are fit than there are aberrations present in the system. In the following, we show experimentally that this is not the case for Method B.
Using the same center position and radius, different amount of Zernike degrees can be fit in order to see the convergence behavior of gathered coefficients. Due to the fact that there are only 1, 2, or 3 Zernike modes present in the specific Zernike experiment, the convergence behavior of these experiments can be visualized.
Based on analysis of the aberrations listed in Table 4.1, the single Zernike experiments 5_1 and 6_4 and 3 random Zernike experiments 3_zerns_1 and 3_zerns_3, one can see in Figure 4.1 and Figure 4.2 that the coefficients seem to be measured independently from each other, and there is no difference between the two methods. However, looking at the subsequent Zernike experiments sub_zerns_1 and sub_zerns_3 a difference can be seen between the two reconstruction methods. The initial guesses of Method A over-or underestimate the presence of the aberration when the fit is performed with a maximum degree of eight. From Figure 4.4 it can be seen that defocus is overestimated by more than 50% until 4 degrees are fit. At maximum degree of four, the spherical aberration coefficient a 0 4 is overestimated. All values seem to be within normal range when fitting a maximum degree of eight. The same can be seen in Figure 4.5, where the coefficient a 3 3 is overestimated at maximum degrees of three and four. Method B does not present this over-estimation.

Conclusion and outlook
In conclusion, we have implemented a new wavefront reconstruction method (Method B) for Shack-Hartmann sensors based on Zernike expansion of derivatives of the Zernike circle polynomials that was introduced in (Janssen, 2014). We have shown with experiments that Method B is advantageous as compared to other known methods such as a least-squares method (Method A) due to the lack of cross-talk between the coefficients. We have also compared Method B to Method A for reconstructing the wavefront using a Shack-Hartmann sensor for some sets of aberrations.
Based on the error values of the reconstruction with optimal amount of Zernike coefficients, the quality of the fit is in general similar for Method A and Method B. When less than the optimal Zernike powers were fit, it was seen that Method B estimates the coefficients more accurately. For the single Zernike aberration case, it was shown that only Method A shows cross-coupling of higher order aberrations, while Method B does not.  As to computational load, we may point out that these are of comparable order, with Method B somewhat more demanding due to the complex arithmetics and the alignment needed for compatibility with the ANSI-format of the Shack-Hartmann sensor.
Finally, we would like to point out that there are also alternative expansions to treat this problem, namely the use of eigenfunctions of the Laplacian with Neumann boundary conditions (Shengyang H and J, 2011), (Shengyang H and J, 2013) that do have orthogonal gradients. It would be interesting as future work to compare this method with Method A and Method B considered in the present paper.

A Relation with orthogonal vector polynomials
In the course of our investigations, we became aware of a seemingly different approach to obtain the wavefront aberration coefficients from wavefront derivative data. We shall describe and relate this approach to our method in the framework and notations of (Janssen, 2014), so that the Zernike circle polynomials are unnormalized and have exponential azimuthal dependence facilitating mathematical developments. The reader will have no particular problems in reformulating the main results of this appendix in terms of the ANSI-style circle polynomials using (Eqs. 2.14−16). Thus, we have for integer n and m such that n − |m| is even and non-negative: with real ν, µ such that ν 2 + µ 2 ≤ 1 and ν + iµ = ρe iθ ; ν = ρ cos θ; µ = ρ sin θ, (A.2) and the radial polynomials R |m| n (ρ) given by Eq. 2.11. We now sketch the approach in (Mahajan and Acosta, 2017), where the notations and conventions differ from our present analysis. The approach (Mahajan and Acosta, 2017) uses the notion of vector polynomials G m n = G m n (ρ, θ) = (G m n,1 (ρ, θ), G m n,2 (ρ, θ)) ∈ C 2 (A.3) for integer n, m such that n − |m| is even and non-negative with n = 0, that satisfy 1 π should be satisfied at the rim ρ 2 = ν 2 + µ 2 = 1 of the pupil. Next, in (Mahajan and Acosta, 2017), subsect. 2.B.1, the G m n are required to be irrotational, meaning that the ∇ × G m n = 0, because of considerations of minimal noise propagation. This condition of irrotationality is satisfied when there is a scalar function U m n on the pupil such that should hold on the rim ρ = 1 of the pupil. In (Mahajan and Acosta, 2017) subsect. 2.B.2-3, the U m n are found by writing the condition Equation A.13 in polar coordinates, with separate consideration of the cases m = 0 and m = 0, and explicitly using the series representation of radial polynomials. This yields the two components of G m n in trigonometric polynomial form, for which it can be shown that the boundary condition Equation A.11 is satisfied as well. In (Janssen, 2014) We shall verify below that, when n = |m| + 2, |m| + 4, ..., this U m n also satisfies the boundary condition Equation A.14, and that a concise formula for G m n = ∇U m n in terms of Zernike circle polynomials results. For the case that n = |m|, we shall show that .14, and we find a concise formula in terms of the Zernike circle polynomials for G m |m| as well. This then can be used to show that α m n of Equation A.9 actually coincide with the LMS estimator found in (Janssen, 2014), sect. 3. The trigonometric/polynomial solution form of G m n found in (Mahajan and Acosta, 2017) involves the coefficients in the series representation of R |m| n (ρ) in Equation A.16, and these become awkward to use when the degree becomes large (n should be limited to ≤ 44 when using double precision). This problem is virtually absent when the representation of the G m n in terms of the Zernike circle polynomials is used since there are nowadays several methods for reliably computing (the radial parts of the) Zernike circle polynomials of arbitrary large degree n and azimuthal order m, (Janssen and Dirksen, 2007), (Shakibaei and Paramesran, 2013).
We now observe that the three series in Equation A.32 have the same terms, except that the second series omits the term l = 0, 1 from the first series and the third series omits the terms with l = 0 from the first series. Using Eq. (A.28), we see that the terms in the first series with l = 2, 3, ...(n+2−|m|)/2 are canceled, and we get In a similar fashion the case n = |m| > 0 can be handled using Equation A.37 and Equation A.38, and this yields that the α m |m| from Equation A.40 and Equation A.41 coincide with the right-hand side integral expression in Equation A.9.