Dwell time for optical fabrication using the modified discrete convolution matrix method

Ximing Liu; Ximing Liu; Ximing Liu; Ximing Liu; Longxiang Li; Longxiang Li; Longxiang Li; Longxiang Li; Xingchang Li; Xingchang Li; Xingchang Li; Xingchang Li; Feng Zhang; Feng Zhang; Feng Zhang; Feng Zhang; Feng Zhang; Xuejun Zhang; Xuejun Zhang; Xuejun Zhang; Xuejun Zhang

doi:10.1364/AO.523682

1. INTRODUCTION

The increasing demand for ultra-precision optical components, such as extreme ultra-violet (EUV) [1], X-ray [2], and lithographic objective [3] mirrors, etc., alongside other precision components, has propelled the rapid advancement of ultra-precision polishing technology. High deterministic polishing techniques, such as magnetorheological finishing [4] and ion beam figuring [5], have been developed to satisfy these demands. Surface removal amount in these techniques is achieved through the convolution of removal tools and dwell time (DT). The analytical method for the DT of a circular mirror is expressed by Eq. (1) as follows:

(1)$$\begin{split} E( x,y )&=R( \omega ,\upsilon )**T( \varepsilon ,\eta ) \\ & = \iint_{{{{\varepsilon }^{2}}+{{\eta }^{2}}\lt R_{t}^{2}}}R( x-\varepsilon ,y-\eta )T( \varepsilon ,\eta ){\rm d}\varepsilon {\rm d}\eta , \end{split}$$

where $E({x,y})$ represents the surface error, $R({\omega ,\upsilon})$ denotes the removal function, $T({\varepsilon ,\eta})$ corresponds to the DT, $**$ represents the convolution, and ${R_t}$ denotes the tool move-area radius.

Many researchers have proposed DT algorithms, broadly categorized into six methods: Bayesian method [5], Fourier transform method [6–10], direct convolution method [11], polynomial fitting method [12–14], genetic algorithms method [15], and discrete convolution matrix method (hereafter referred to as the original method). The Bayesian method ensures non-negativity; however, it is limited to circular symmetry removal functions [16]. The Fourier transform method offers rapid solutions but necessitates a threshold factor, rendering the solution unstable. The direct convolution method, while flexible, cannot guarantee non-negativity of the DT [16]. The polynomial fitting method ensures smoothness; however, it relies on the choice of polynomials for accuracy. The genetic algorithm method relatively consumes much calculation time [17].

Finally, the original method, employing a matrix approach, is flexible and insensitive to path clutter degree. The tool path is typically discretized for computational purposes. Commonly employed discretization paths include grating paths [18], spiral paths [19], Peano paths [20], and other random paths [21–25]. Because the random path can produce less deterministic structure on the surface, it has been widely studied. Thus, it is widely used in a discrete form [17,26–32], as expressed in Eqs. (2)–(4):

Fig. 1. (a) Uniform isometric discrete grid. (b) Grating path induced by (a). (c) Hilbert path induced by (a).

Download Full Size | PDF

(2)$$E({{x_i},{y_i}} ) = \sum\limits_j {R({{x_i} - {\varepsilon _j},{y_i} - {\eta _j}} )T({{\varepsilon _j},{\eta _j}} )} ,$$

(3)$${e_i} = E({{x_i},{y_i}} ),{r_{\textit{ij}}} = R({{x_i} - {\varepsilon _j},{y_i} - {\eta _j}} ),{t_j} = T({{\varepsilon _j},{\eta _j}} ),$$

(4)$$\left[{\begin{array}{*{20}{c}}{{r_{11}}}&{{r_{12}}}& \cdots &{{r_{1N}}}\\{{r_{21}}}&{{r_{22}}}& \cdots &{{r_{2N}}}\\ \vdots & \vdots & \ddots & \vdots \\{{r_{M1}}}&{{r_{M2}}}& \cdots &{{r_{\textit{MN}}}}\end{array}} \right]\left[{\begin{array}{*{20}{c}}{{t_1}}\\{{t_2}}\\ \vdots \\{{t_N}}\end{array}} \right] = \left[{\begin{array}{*{20}{c}}{{e_1}}\\{{e_2}}\\ \vdots \\{{e_M}}\end{array}} \right] \Leftrightarrow {R^\prime _o}{\textbf{T}^\prime_o}= \textbf{E}^\prime,{\rm s.t.}{T^\prime _o} \ge 0,$$

where ${t_j} = T({{\varepsilon _j},{\eta _j}})$, $1 \le j \le N$ represents the DT of vectorized path points; $N$ represents the number of path points; ${e_i} = E({{x_i},{y_i}})$; $1 \le i \le M$ represents the surface points of surface error; $M$ represents the number of surface points; and ${r_{\textit{ij}}} = R({{x_i} - {\varepsilon _j},{y_i} - {\eta _j}})$ represents the $({i,j})$ matrix element of the convolution kernel. They are simply called as vectors ${R^\prime _o} = [{{r_{\textit{ij}}}}],\textbf{E}^\prime = [{{e_i}}],{\textbf{T}^\prime_o} = [{{t_j}}]$$,1 \le i \le M,1 \le j \le N$.

Further, calculated surface residuals are Eq. (5) referred to [33]

(5)$${\boldsymbol {\Delta E}} = \left({\textbf{E}^\prime - {{R^\prime}_o}{{{\textbf{T}^\prime}}_o}} \right).$$

Nonetheless, the analytic expression Eq. (1) of the DT problem, coupled with the uniform isometric discrete grid discretization form, is elucidated as in Eq. (6). The uniform isometric discrete form is illustrated in Fig. 1(a). The grating path [Fig. 1(b)] and the Hilbert path [Fig. 1(c)] are the outcomes of the same discretization. Until their discrete grids are consistent, the result of calculating the DT using Eq. (4) is consistent, provided the grid is not sufficiently dense to simulate the tool marks error [34].

However, depending on the original method can cause problems that the absence of the $\Delta \varepsilon \Delta \eta$ term in Eq. (2). In practice, we found that the term can reduce the uniformity of surface residuals. So, we came up with Eq. (6) to clarify the effect of the term:

(6)$$\begin{split}E({{x_i},{y_i}} ) &\approx \mathop {\lim}\limits_{\stackrel{M \to \infty }{N \to \infty }} \sum\limits_{m = 1}^{{N_\varepsilon}} {\sum\limits_{n = 1}^{{N_\eta}} {R({{x_i} - {\varepsilon _m},{y_i} - {\eta _n}} )T({{\varepsilon _m},{\eta _n}} )}} \Delta \varepsilon \Delta \eta \\& \approx \sum\limits_l^{{N_\varepsilon}{N_\eta}} {R({{x_i} - {\varepsilon _l},{y_i} - {\eta _l}} )T({{\varepsilon _l},{\eta _l}} )} \Delta {\varepsilon _l}\Delta {\eta _l},\end{split}$$

where ${N_\varepsilon}$ represents the number of discretization points along $\varepsilon$, $N_\eta$ represents the number of discretization points along $\eta$, and $l$ denotes the number of the dwell points in path order.

This study principally addressed these terms via analysis and simulation. The remainder of this paper is structured as follows: In Section 2, the modified method of the discrete convolution matrix method is introduced. In Section 3, a detailed analysis spiral path weight is presented and simulation results involving surface residuals of both the spiral path and grating-like path are displayed. The conclusions are presented in Section 4.

2. MODIFIED DISCRETE CONVOLUTION MATRIX METHOD

When the physical meaning of weight is the area represented by discrete points, an approximate general formula can be obtained.

Consequently, Eq. (6) is rewritten as Eq. (7), which is illustrated to show in Fig. 2(a):

(7)$$\begin{split}E({{x_i},{y_i}} ) &= \sum\limits_l {R({{x_i} - {\varepsilon _l},{y_i} - {\eta _l}} )T({{\varepsilon _l},{\eta _l}} )} {S_l}\\ &\Rightarrow \textbf{E}^\prime = {{R}_m^\prime}{{\textbf{T}}_m^\prime},\end{split}$$

where ${S_l}$ represents the partitioning area (or weight of $({{\varepsilon _l},{\eta _l}})$).

Fig. 2. Schematic diagram of the contribution value of the $l$-th dwell point to the $i$-th removal amount. (a) One of the weights of VP partitioning of domain in 3D. (b) Top view of (a).

Download Full Size | PDF

The modified matrix element should be as in Eq. (8):

(8)$${{r}_{\textit{il}}^\prime} = R\left({{x_i} - {\varepsilon _l},{y_i} - {\eta _l}} \right){S_l},\quad{{R}_m^\prime} = \left[{{{r}_{\textit{il}}^\prime}} \right].$$

The theoretical modification of the Cartesian coordinate system is as expressed in Eq. (9), wherein the term $\Delta \varepsilon \Delta \eta$ is added:

(9)$${S_l} = \Delta {\varepsilon _l}\Delta {\eta _l}.$$

The theoretical modification of the polar coordinate system is presented in Eq. (10), wherein the term $\rho \Delta \rho \Delta \theta$ is added [35]:

(10)$${S_l} = {\rho _l}\Delta {\rho _l}\Delta {\theta _l}.$$

To calculate the general non-uniform discrete grid, we considered the area represented by it as the weight. The area of a Voronoi polygon (VP) [36] was considered as its weight via VP partitioning.

Fig. 3. (a) Random discretized grid. (b) Voronoi diagram according to (a). (c) Grating-like path induced by (a).

Download Full Size | PDF

The distance of any point in ${\textit{VP}_l}$ to the control point $({{\varepsilon _l},{\eta _l}})$ is less than the distance to $({{\varepsilon _k},{\eta _k}}),k \ne l$. Here, $\textit{VP}_l$ is the $l $-th Voronoi equivalence class, and $({{\varepsilon _l},{\eta _l}})$ is the representative of ${\textit{VP}_l}$. Random points are generated to compute the DT problem, as shown in Fig. 3(a). The VP is computed as in Fig. 3(a) using the MPT toolbox [37], which is shown in Fig. 3(b). Machine dynamic performance was considered to render the path as straight as possible, as shown in Fig. 3(c). The area of ${\textit{VP}_l}$ it represents was calculated as Eq. (11):

(11)$${S_l} = {{\rm Area}} \left({{\textit{VP}_l}} \right),$$

where ${\rm Area} (\cdot)$ is a function of the area of a domain. An example is shown in Fig. 2(b).

Remark: The discrete points of partitioning of Eqs. (9) and (10) are the vertices of the partitioned geometries and an example is shown Fig. 4(a); however, in the partitioning of the VP, they are representatives of the Voronoi equivalence classes, and an example is shown in Fig. 4(b).

Fig. 4. Two kinds of weights for the same grid. (a) Polar differential area weight referred to Eq. (10). (b) VP area weight referred to Eq. (11).

Download Full Size | PDF

The DT calculated according to Eq. (7) cannot be used for CNC code writing directly. There is another step such as Eq. (12) weight allocation. Because weights are added to the convolution kernel matrix ${R^\prime _m}$ during the calculation, but the actual removal function matrix $R({\omega ,\upsilon})$ does not obtain additional weights. The calculated weights are reassigned to the DT to match the actual removal function for correct and accurate material removal:

(12)$$\begin{split}{\textbf{T}_{{\rm actual}}} = {{\textbf{T}^\prime_m}} \times {\rm diag}(S ),\\{\rm where}\quad S = \left[{{S_1},{S_1}, \ldots ,{S_N}} \right],\\{\rm diag}(S ) = \left[{\begin{array}{*{20}{c}}{{S_1}}&0&0&0\\0&{{S_2}}&0&0\\0&0& \ddots &0\\0&0&0&{{S_N}}\end{array}} \right].\end{split}$$

Thus, the residual is calculated as Eq. (13):

(13)$${\boldsymbol {\Delta E}} = \left({\textbf{E}^\prime - {{R}_o^\prime}{{\textbf{T}}_{\rm{actual}}^\prime}} \right).$$

The flow chart of the algorithm is given, as shown in Fig. 5. $\textit{mpt}\_\textit{voronoi}(*)$ represents the function for computing the VP with boundaries, derived from MPT toolbox. $bd$ denotes the artificially set Voronoi boundary; the boundary used in this paper is defined as the optical component aperture. $\textit{diag}(*)$ means to make a diagonal matrix with the parameter vector. $\textit{CNC}(*)$ stands for writing computer numerical control codes for machine tools based on the parameter dwell time vector. $\textit{mef}$ represents machine executable file.

Fig. 5. Modified dwell time algorithm flow chart.

Download Full Size | PDF

3. ANALYTIC AND NUMERICAL EXAMPLE

A. Spiral Path

First, the spiral path (emphasizing that the direction of the removal function remained constant with the dwell point) example was calculated, and its discrete points distribution is shown in Fig. 6(a). The VP diagram of the spiral path is also shown in Fig. 6(b), and the continuous path is illustrated in Fig. 6(c).

Fig. 6. (a) Spiral discrete grid. (b) Voronoi diagram according to (a). (c) Spiral path induced by (a).

Download Full Size | PDF

The spiral path expression is expressed as Eq. (14):

(14)$$\rho (\theta ) = a + \frac{b}{{2\pi}}\theta ,\quad\theta \in [{0,2\pi C} ],$$

where $a$ denotes the inner diameter of the mirror. When the mirror to be machined lacked a central hole, $a = 0$. Further, $b$ denotes the pitch of the spiral path, $C = {{\rm ceil}} [({{R_e} - a})/b]$ represents the maximum integer turns of the spiral path, ${\rm ceil}(\cdot)$ and corresponds to the ceiling function.

This integral problem should be computed using Eq. (7) in the form of Eq. (10). The number of points in each circle should vary to ensure uniformity among discrete points. We set the $q $-th circle to be dispersed by ${\gamma _q}$ points [Fig. 6(c)]. Further, we set ${\gamma _0} = 0$. First, the removal amount for any surface error point on the $q $-th circle was calculated:

(15)$$\begin{split}{E_q}({x,y} ) &= {({2\pi} )_q}\sum\limits_{{n_q} = 1}^{{\gamma _q}} {R\left({x - {\varepsilon _{{n_q}}},y - {\eta _{{n_q}}}} \right)} \\&\quad\times T\left({{\varepsilon _{{n_q}}},{\eta _{{n_q}}}} \right) \times\left({a + \frac{b}{{2\pi}}{\theta _{{n_q}}}} \right) \times \frac{b}{{2\pi}}\Delta {\theta _q}\\{\rm where},\quad {\theta _{{n_q}}} &= \frac{{2\pi ({{n_q} - 1} )}}{{{\gamma _q}}} + 2\pi ({q - 1} ),\\\Delta {\theta _q} &= \frac{{2\pi}}{{{\gamma _q}}},1 \le {n_q} \le {\gamma _q},\end{split}$$

where $({{\varepsilon _{{n_q}}},{\eta _{{n_q}}}})$ denotes the ${n_q}$-th dwell point of the $q $-th circle. Considering the interdependence between $\rho$ and $\theta$, discretization was applied solely to $\theta$ and acts on $\rho$ as $d\rho \to \Delta {\rho _q} = {{b\Delta {\theta _q}} / {({2\pi})}}$. ${({2\pi})_q} = {\int_{2\pi ({q - 1})}^{2\pi q}} {\rm d}\theta = \sum\nolimits_{{n_q} = 1}^{{\gamma _q}} \Delta {\theta _q}$ represents the integral contribution of the polar angle $\theta$ of the $q $-th circle.

Consequently, the first multiple integral of $\theta$ served as an approximation to $2\pi$. The discretization of $\theta$ was only applicable to the second multiple integral of $\rho$. The source of $2\pi$, denoted as ${({2\pi})_q}$, must be identified. Moreover, the approximation error under condition $b \le {R_e}/50$ was disregarded by empiricism, where ${R_e}$ represents the radius of the component.

The total removal amount calculated by all the dwell points is the summation of the contribution from each circle:

(16)$$\begin{split}E({x,y} ) & = \sum\limits_{q = 1}^C {{E_q}({x,y} )} \\ & = \sum\limits_{q = 1}^C {{({2\pi} )}_q}\sum\limits_{{n_q} = 1}^{{\gamma _q}} {R\left({x - {\varepsilon _{{n_q}}},y - {\eta _{{n_q}}}} \right)} \\&\quad\times T\left({{\varepsilon _{{n_q}}},{\eta _{{n_q}}}} \right) \times \left({a + \frac{b}{{2\pi}}{\theta _{{n_q}}}} \right) \times \frac{b}{{2\pi}}\Delta {\theta _q} \\ & = \sum\limits_{l = 1}^N {R\left({x - {\varepsilon _l},y - {\eta _l}} \right)} \times T\left({{\varepsilon _l},{\eta _l}} \right) \times {({2\pi} )_q} \\&\quad\times \left({a + \frac{b}{{2\pi}}{\theta _l}} \right) \times \frac{b}{{2\pi}}\Delta {\theta _q}\\ {\rm where} \quad 1 \le l & = {n_q} + \sum\limits_{\alpha = 0}^{q - 1} {{\gamma _\alpha}} \le N = \sum\limits_{q = 1}^C {\sum\limits_{{n_q} = 1}^{{\gamma _q}} 1} .\end{split}$$

Consequently, the matrix element weight can be represented as Eq. (17):

(17)$$\begin{split}{S_l} & = {({2\pi} )_q}\left[{a + b \times \frac{{({{n_q} - 1} )}}{{{\gamma _q}}} + b \times ({q - 1} )} \right]\\\frac{b}{{{\gamma _q}}}& = \frac{{{\alpha _1}}}{{{\gamma _q}}} + \frac{{{\alpha _2}({{n_q} - 1} )}}{{\gamma _q^2}},\\ {\rm where}\quad {\alpha _1} & = {({2\pi} )_q}\left[{ab + {b^2} \times ({q - 1} )} \right],\quad {\alpha _2} = {({2\pi} )_q}{b^2}.\end{split}$$

Numerical calculations are used to build intuition. The radius of the component is ${R_e} = 50\;{\rm mm}$ and the surface error employed in the simulation is depicted in Fig. 7(a). To prevent edge effects, a clear aperture (CA) was set for the surface residuals observation. The radius of the CA is ${R_{\rm{CA}}} = 40\;{\rm mm}$. The surface error was sampled at a rate of 1 mm per point using a uniform isometric discrete grid. The removal function used is shown in Fig. 7(b). The expression for the spiral path was $\rho (\theta) = [{{1 / {({2\pi})}}}] \times \theta ,\theta \in [{0,2\pi \times 60}]$. The path coverage ${R_t} = C \times b = 60\;{\rm mm}$ is larger than the component diameter to suppress the ringing effect from deconvolution [5]. The initial discrete parameter was ${\gamma _q} = 10 \times q$. According to these parameters, Eq. (17) can be written as (18):

Fig. 7. (a) Surface error. (b) Removal function. (c) DT calculated using Eq. (4). (d) DT calculated using Eq. (12). (e) The residuals calculated using Eq. (5) with DT (c). (f) The residuals calculated using Eq. (13) with DT (d).

Download Full Size | PDF

(18)$$\begin{split}{\alpha _1} & = {({2\pi} )_q} \times ({q - 1} ),\quad {\alpha _2} = {({2\pi} )_q},\\ l & = {n_q} + \sum\limits_{\alpha = 0}^{q - 1} {10 \times \alpha} = {n_q} + 5 \times q \times ({q - 1} ),\\{S_l} & = \frac{{{{({2\pi} )}_q} \times ({q - 1} )}}{{10q}} + \frac{{{{({2\pi} )}_q} \times \left({{n_q} - 1} \right)}}{{100{q^2}}},\quad 1 \le {n_q} \le 10 \times q.\end{split}$$

Figure 7(c) shows the DT calculated using the original method Eq. (4), and Fig. 7(d) presents the DT calculated using the modified method Eq. (12). The results of VP weights calculated under the above conditions are approximately identical with Fig. 7(d). The oval ring structure and radiation stripes on the DT in Figs. 7(c) and 7(d) resulted from the unavoidable ringing effect and interaction of uniform isometric on spiral grids. There are values outside the aperture because the calculation is extended, and these small removal amounts also determine the residual amount of the surface in the actual processing, which will load a retaining edge fixture outside the component, so that the removal function continues to remove this area.

From the color bar, it can be seen that compared with the original method, the dwell time value of most dwell points in the modified method are still relatively distributed basically the same. The maximum dwell time of the modified method is significantly smaller than that of the original method. The obvious difference is that the dwell time center calculated by the modified method has a dark spot compared with the original method.

Figures 7(e) and 7(f) show the residuals calculated using the original and modified polar weight methods. The ring structure is inherited from the dwell time. In the residuals calculated by the original method, a removal pit appears obviously. The uniformity parameters peak to valley (PV) and root-mean-square (RMS) values of the calculated residuals only approach $\lambda /20$ and $\lambda /145$. However, by the modified method, the PV and RMS values of the calculated residuals only approach $\lambda /40$ and $\lambda /190$, which has been relatively improved by 51.4% and 23.1%. It shows that the modified method improves the uniformity and accuracy of residuals in component polishing. The convergence rate was $\eta = ({\rm RMS}_{\rm error}-{\rm RMS}_{\rm residual})/{\rm RMS}_{\rm error}\in [0,1]$, which relat-ively optimized by 0.78%.

Subsequently, the mechanism of the modified method is further analyzed. The three weight curves, as shown in Fig. 8(a), were used to illustrate the cause of the phenomenon in Figs. 7(d) and 7(f). The first observation is that the curves calculated by Eqs. (10) and (11) are almost the same, because they are the two representations of the area weight. The weight of the original method is set to 1 for normalization purposes, and it can be set to any value desired to represent that they have equal power to participate in the calculation. However, if the equal weight is not 1, the $\textit{diag}(*)$ transformation by Eq. (12) is required.

Fig. 8. (a) Weights for three different methods. (b) Local magnification of purple box.

Download Full Size | PDF

It can be seen from Eq. (18) that the weight is a linear increasing region for each fixed ${\gamma _q}$ that is equivalent to the $q $-th circle. The slope of the linear domain is proportional to ${1 / {\gamma _q^2}}$, and since the iteration rule increase ${\gamma _q}$, the slope of each linear region decreases with $\rho$ increased. It eventually converges to a constant value ${{2\pi} / {10}}$, which is given by Eq. (19) referred to as the Sandwich theorem when $q$ increased as equivalent as $\rho$ increased. The upper and lower limits $({ul \;\& \;ll})$ that vary with $q$ are its envelopes. The eventual result of a fixed value means that the weights of the peripheral points are almost equal, which is why most of the DT distribution of dwell points did not differ. And since the final equal weight is ${{2\pi} / {10}} \lt 1$, according to the above $\textit{diag}(*)$ transformation, the maximum value of the improved method is less than the original method:

(19)$$\begin{split}ll & = \frac{{{{({2\pi} )}_q} \times ({q - 1} )}}{{10q}} + 0 \le {S_l} \le ul\\&= \frac{{{{({2\pi} )}_q} \times ({q - 1} )}}{{10q}} + \frac{{{{({2\pi} )}_q} \times ({10q - 1} )}}{{100{q^2}}},\\ \mathop {\lim}\limits_{q \to \infty} ({ll} )& = \mathop {\lim}\limits_{q \to \infty} ({ul} ) = \frac{{{{({2\pi} )}_q}}}{{10}}.\end{split}$$

However, the calculations vary dramatically near the center, which is magnified in Fig. 8(b). Compared with the original method, the improved method reduces the enthusiasm of the near central point to participate in the calculation. Observing Fig. 6(a), it can be seen that the path points in the center are very dense, so the area they can represent is very small, and the center point can only represent itself in the calculation (i.e., weight is 0), which is the reason why the DT center calculated by the modified method is a dark spot. It will prevent the values at these points from being innocently apportioned, making the DT calculation more accurate. It is not difficult to analyze the modified residuals due to the small DT (i.e., small removal) at the center point, resulting in uniform residuals without a deep removal pit in the center, compared with the original method. This holds significance in optical component fabrication that necessitates the high uniform convergence of error.

B. Grating-Like Path

First, the grating-like path (emphasizing that the direction of the removal function remained constant with the dwell point) example is shown in Fig. 3(c). The radius of the component is ${R_e} = 50\;{\rm mm} $ and surface error employed in the simulation is depicted in Fig. 9(a). The removal function used is shown in Fig. 9(b).

Fig. 9. Simulation data of grating-like path. (a) Surface error. (b) Removal function. (c) DT calculated using Eq. (4). (d) DT calculated using Eq. (12). (e) The residuals calculated using Eq. (5) with DT (c). (f) The residuals calculated using Eq. (13) with DT (d).

Download Full Size | PDF

Figures 9(c) and 9(d) show the DT calculated using the original method and the modified method with VP weight form referred to Fig. 3(b) and Eq. (11). Compared with the original method, the total DT of the modified method is much the same, but the PV value is twice that of the original method. In order to explain the DT distribution, the discrete path points are also scattered in the diagram. Both methods show that the DT in the region with low dwell points density is higher than that in the surrounding region, which is more prominent by the modified method because the VP area (i.e., the weight) represented at this dwell point is very large so that more DT is allocated here under the almost same total DT. In the original method, each point participates in the calculation with the same weight such that it pulls down the contribution of dwell point at low dwell points density, so the DT is smaller.

Figures 9(e) and 9(f) show the residuals calculated using the original and modified VP weight methods. The PV and RMS values have been relatively deteriorated by 1.2% and 16.5%. The convergence rate $\eta$ relatively deteriorated by 0.5%. Compared with the original method, the modified method reduces the convergence rate and makes the PV and RMS of the surface decrease, because under the same surface error and removal function, the DT distribution with large amplitude fluctuation means that the removal amount distribution is more non-uniform, resulting in more chaotic distribution of residuals. However, we believe that DT calculation by the modified method is the theoretical upper limit of the PV, RMS, and convergence rate obtained by this group of the parameters of surface error, removal function, and discrete dwell points participating in the DT calculation. Meanwhile, the specific numerical values are also depended on the specific optimization programming algorithm. It is unrealistic about the high convergence rate calculated by the original method, due to the above analysis of the DT calculation. There are many error sources in polishing tasks, and this tiny method error is usually neglected in previous tasks. The proposed modified method can help researchers more confidently conduct research to find other sources of error in polishing.

In short, the weight of the modified method is derived from the area represented by the dwell point, which is sensitive to the distribution density of the dwell points. Therefore, for the tool path with the non-uniform discretization, it is necessary to consider the modified method to achieve the maximum accuracy of the DT calculation.

4. CONCLUSIONS

This study proposed a weight modified matrix element method, and the modification of the deconvolution matrix element was obtained. A detailed analytical calculation and simulation of the spiral path correction matrix elements were conducted. The simulation results indicated that the modified method improved the uniformity of the calculation solution of the spiral path of dwell time, and the weight curve was used to analyze how the modified method rendered the uniform. The necessity of this modified method for accurate dwell time calculation under non-uniform discretization was further substantiated through grating-like path simulation. The study offers valuable guidance for accurate dwell time calculation in achieving uniform error converges during the polishing process.

Although this modification has notably improved the error convergence and accuracy in solving dwell time for non-uniform discretization, we acknowledge that key approximations with specific conditions still exist. Consequently, the application of this calculation method may face limitations when encountering a large sparse path (i.e., the dwell points density is very low), such as the calculation of scenarios where the approximation of ${({2\pi})_q}$ remains unsatisfied. In addition, Voronoi polygon calculation of a large number of discrete point data is relatively slow, although the method improves the accuracy, it is bound to reduce some efficiency. In the future, we will focus on improving the Voronoi polygon generation algorithm to achieve fast and high-precision dwell time calculation of the non-uniform discrete tool path.

Funding

National Natural Science Foundation of China (62275246; 62075218; 12003034); Youth Innovation Promotion Association of the Chinese Academy of Sciences (2021215); National Key Research and Development Program of China (2022YFB3403405).

Acknowledgment

The authors would like to thank the anonymous referees for their valuable suggestions and comments.

Disclosures

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

REFERENCES

1. L.-H. Wang, X.-K. Wang, and B. Chen, “Study for dual-function EUV multilayer mirror,” Opt. Laser Technol. 40, 571–574 (2008). [CrossRef]

2. S. Ichimaru, M. Hatayama, T. Ohchi, et al., “A proposal for a collecting mirror assembly for large divergence x-ray sources,” Appl. Opt. 53, 7348–7354 (2014). [CrossRef]

3. J. Sun, X. Li, W. Zhang, et al., “High-reflectivity mirrors by Al₂O₃, LaF₃ and AlF₃ for 193nm application,” Opt. Laser Technol. 56, 65–70 (2014). [CrossRef]

4. B. Wang, F. Shi, W. Zhang, et al., “Key technology research on magnetorheological finishing based on suppression of surface mid-spatial frequency ripple errors,” Opt. Mater. Express 12, 3213–3224 (2022). [CrossRef]

5. C. Jiao, S. Li, and X. Xie, “Algorithm for ion beam figuring of low-gradient mirrors,” Appl. Opt. 48, 4090–4096 (2009). [CrossRef]

6. S. Wan, X. Zhang, M. Xu, et al., “Region-adaptive path planning for precision optical polishing with industrial robots,” Opt. Express 26, 23782–23795 (2018). [CrossRef]

7. T. Wang, L. Huang, H. Kang, et al., “RIFTA: A Robust Iterative Fourier Transform-based dwell time algorithm for ultra-precision ion beam figuring of synchrotron mirrors,” Sci. Rep. 10, 8135 (2020). [CrossRef]

8. T. Wang, L. Huang, H. Choi, et al., “RISE: robust iterative surface extension for sub-nanometer X-ray mirror fabrication,” Opt. Express 29, 15114–15132 (2021). [CrossRef]

9. T. Wang, L. Huang, M. Vescovi, et al., “Universal dwell time optimization for deterministic optics fabrication,” Opt. Express 29, 38737–38757 (2021). [CrossRef]

10. S. Tanebe and S. I. Satake, “Iterative blind deconvolution method for dwell-time adjustment,” Appl. Opt. 61, 1743–1749 (2022). [CrossRef]

11. C. Wang, W. Yang, Z. Wang, et al., “Dwell-time algorithm for polishing large optics,” Appl. Opt. 53, 4752–4760 (2014). [CrossRef]

12. W. Zhu and A. Beaucamp, “Zernike mapping of optimum dwell time in deterministic fabrication of freeform optics,” Opt. Express 27, 28692–28706 (2019). [CrossRef]

13. Y. Wang, Y. Zhang, R. Kang, et al., “An elementary approximation of dwell time algorithm for ultra-precision computer-controlled optical surfacing,” Micromachines 12, 471 (2021). [CrossRef]

14. Z. Li, R. Wang, X. Zhang, et al., “B-spline surface approximation method for achieving optimum dwell time in deterministic polishing,” J. Mater. Process. Technol. 318, 118031 (2023). [CrossRef]

15. J. González-García, A. Cordero-Dávila, I. Leal-Cabrera, et al., “Calculating petal tools by using genetic algorithms,” Appl. Opt. 45, 6126–6136 (2006). [CrossRef]

16. M. Xie, J. Cheng, Z. Ma, et al., “Review of dwell-time algorithms for magnetorheological finishing,” J. Adv. Manuf. Technol. 130, 3151–3178 (2024). [CrossRef]

17. A. Santiago-Alvarado, J. González-García, C. Castañeda-Roldan, et al., “Use of linear programming to calculate dwell times for the design of petal tools,” Appl. Opt. 46, 4642–4649 (2007). [CrossRef]

18. M. Rososhansky and F. Xi, “Coverage based tool-path planning for automated polishing using contact mechanics theory,” J. Manuf. Syst. 30, 144–153 (2011). [CrossRef]

19. V. S. Negi, T. Wang, H. Garg, et al., “Random adaptive tool path for zonal optics fabrication,” Opt. Express 30, 29295–29309 (2022). [CrossRef]

20. H. Y. Tam, H. Cheng, and Z. Dong, “Peano-like paths for subaperture polishing of optical aspherical surfaces,” Appl. Opt. 52, 3624–3636 (2013). [CrossRef]

21. K. Takizawa and A. Beaucamp, “Comparison of tool feed influence in CNC polishing between a novel circular-random path and other pseudo-random paths,” Opt. Express 25, 22411–22424 (2017). [CrossRef]

22. H. Li, X. Li, S. Wan, et al., “High-efficiency smooth pseudo-random path planning for restraining the path ripple of robotic polishing,” Appl. Opt. 60, 7732–7739 (2021). [CrossRef]

23. C. Wang, Y. Han, H. Zhang, et al., “Suppression of mid-spatial-frequency waviness by a universal random tree-shaped path in robotic bonnet polishing,” Opt. Express 30, 29216–29233 (2022). [CrossRef]

24. T. Wang, H. Cheng, W. Zhang, et al., “Restraint of path effect on optical surface in magnetorheological jet polishing,” Appl. Opt. 55, 935–942 (2016). [CrossRef]

25. C. Wang, Z. Wang, and Q. Xu, “Unicursal random maze tool path for computer-controlled optical surfacing,” Appl. Opt. 54, 10128–10136 (2015). [CrossRef]

26. L. Li, D. Xue, W. Deng, et al., “Positive dwell time algorithm with minimum equal extra material removal in deterministic optical surfacing technology,” Appl. Opt. 56, 9098–9104 (2017). [CrossRef]

27. L. Zhang, Q. Zhao, and C. Fan, “Dwell time algorithm in deterministic polishing of a free-form surface based on the continuous tool influence function,” Appl. Opt. 60, 2704–2715 (2021). [CrossRef]

28. H. Kang, T. Wang, H. Choi, et al., “Genetic algorithm-powered non-sequential dwell time optimization for large optics fabrication,” Opt. Express 30, 16442–16458 (2022). [CrossRef]

29. L. Li, W. Deng, B. Zhang, et al., “Dwell time algorithm for large aperture optical element in magnetorheological finishing,” Acta Opt. Sin. 34, 0522001 (2014). [CrossRef]

30. C. Lv, L. Zou, Y. Huang, et al., “A novel toolpath for robotic adaptive grinding of extremely thin blade edge based on dwell time method,” IEEE ASME Trans. Mechatron. 27, 4429–4439 (2022). [CrossRef]

31. A. Chernyshev, N. Chkhalo, I. Malyshev, et al., “Matrix based algorithm for ion-beam figuring of optical elements,” Precis. Eng. 69, 29–35 (2021). [CrossRef]

32. T. Huang, D. Zhao, and Z.-C. Cao, “Trajectory planning of optical polishing based on optimized implementation of dwell time,” Precis. Eng. 62, 223–231 (2020). [CrossRef]

33. X. Liu, L. Li, X. Li, et al., “Tool mark prediction on the surface of large-aperture mirrors via magnetorheological finishing,” Opt. Express 32, 11150–11170 (2024). [CrossRef]

34. T. Wang, X. Ke, L. Huang, et al., “Computer-controlled finishing via dynamically constraint position-velocity-time scheduler,” J. Manuf. Processes 87, 97–105 (2023). [CrossRef]

35. F. Shi, Y. Dai, X. Peng, et al., “Dwell time algorithm based on vector for MRF process of optics,” J. Natl. Univ. Def. Technol. 31, 103–106 (2009).

36. S. Wan, X. Zhang, W. Wang, et al., “Edge control in precision robotic polishing based on space-variant deconvolution,” Precis. Eng. 55, 110–118 (2019). [CrossRef]

37. M. Kvasnica, P. Grieder, and M. Baoti, “Multi-parametric toolbox (MPT),” 2004,http://control.ee.ethz.ch/∼mpt/.

Dwell time for optical fabrication using the modified discrete convolution matrix method

Abstract

1. INTRODUCTION

2. MODIFIED DISCRETE CONVOLUTION MATRIX METHOD

3. ANALYTIC AND NUMERICAL EXAMPLE

A. Spiral Path

B. Grating-Like Path

4. CONCLUSIONS

Funding

Acknowledgment

Disclosures

Data availability

REFERENCES

Data availability

Cited By

Figures (9)

Equations (19)

Applied Optics