The present group previously proposed a method for determining the first-order derivative matrix (i.e., the Jacobian matrix) of a skew ray by taking all the independent variables of the optical system as the system variable vector, ${\bar{{\textrm{X}}}_{\textrm{sys}}}$. However, many trigonometric function calls, divisions, multiplications and additions were required to determine the ray Jacobian matrix with respect to ${\bar{{\textrm{X}}}_{\textrm{sys}}}$. Accordingly, in the present study, the angular variables in the system variable vector, ${\bar{{\textrm{X}}}_{\textrm{sys}}}$, are replaced with their respective cosine and sine trigonometric functions. The boundary variable vector, ${\bar{{\textrm{X}}}_{\textrm{i}}}$, is similarly redefined such that it includes no angular variables. The proposed method has three main advantages over that previously reported: 1) it is valid for any pose matrix, irrespective of the order in which the rotation and translation motions of a boundary are assigned; 2) it involves only polynomial differentiation, and is thus easily implemented in computer code; and 3) the computation speed of ${{\partial {{\bar{{\textrm{X}}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{{\textrm{X}}}}_{\textrm{i}}}} {\partial {{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}$ is improved by a factor of approximately ten times.
1. Introduction
When optimizing and analyzing optical systems, it is necessary to calculate the Jacobian matrices of various optical quantities (e.g., the exit ray, the optical path difference, and so on) such that their effects on the merit function of the system can be evaluated. In performing such analyses, the relative effect of each variable is determined by computing the partial derivatives of the merit functions of interest with respect to the variable in question. This process can be performed using either Finite Difference (FD) approximation methods [1] or differential methods [2–12]. To the best of the authors’ knowledge, commercial optical analysis and design software packages (e.g., ZEMAX [13] and Code-V [14]) all adopt the former approach to estimate the Jacobian matrix. In the FD approach, each variable is adjusted incrementally whilst the remaining variables are held constant, and the corresponding effects on the optical quantities of interest are evaluated via a raytracing approach. The FD method is relatively straightforward and is easily implemented in computer code. However, the accuracy of the results is reliant on a suitable choice of the step size used to adjust the variables during the tuning stage. For example, an overly-large step size violates the assumption of local linearity, while an overly-small step size reduces the difference between the original solution and the perturbed solution and hence leads to higher rounding errors.
Although several anonymous referees have claimed the ability to compute the Jacobian matrix of any desired merit function with acceptable accuracy using FD methods, most researchers have persisted in the use of differential methods [2–12]. However, the methods presented in [2–12] do not take account of all the independent system variables (e.g., the six pose variables of an element in 3D space). Accordingly, in a previous study by the current group [15], a method was proposed for determining the Jacobian matrix of a skew ray ${\bar{\textrm {R}}}_{\textrm{i}}$ by taking all the independent variables as the system variable vector, ${\bar{{\textrm{X}}}_{\textrm{sys}}}$. However, the method in [15] is valid only for pose matrices specified by translational motion (Eq. (34d) of Appendix A) followed by roll-pitch-yaw (RPY) motion (Eq. (34e) of Appendix A). Thus, many difficulties are encountered when attempting to apply the proposed methodology to pose matrices specified with other combinations of the motion matrices (e.g., RPY motion followed by translation motion). In particular, it is found that the computation of ${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \partial }} \right.} \partial }{\bar{{\textrm{X}}}_{\textrm{sys}}}$, where ${}^0{\bar{{\textrm {A}}}_{\textrm{i}}}$ is the pose matrix of the ith boundary surface relative to the world coordinate frame, involves many trigonometric function calls, divisions, multiplications and additions, as shown in the appendices of Chapter 8 of [15]. In order to overcome these difficulties, the present study replaces the angular variables in the system variable vector, ${\bar{{\textrm{X}}}_{\textrm{sys}}}$, with their respective cosine and sine trigonometric functions. It is shown that the proposed method is applicable to any pose matrix specified by an arbitrary combination of translation and rotation motions. Furthermore, its computation involves only polynomial differentiation, i.e., trigonometric function calls are not required. As a result, the computation speed of ${{\textrm {d}{{\bar{{\textrm{X}}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{d{{\bar{{\textrm{X}}}}_{\textrm{i}}}} {d{{\bar{{\textrm{X}}}}_{sys}}}}} \right.} {\textrm{d}{{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}$ is around ten times faster than that of the method in [15]. The proposed methodology can be extended to higher order ray derivatives to investigate the ray aberrations [16].
2. System variable vector
Consider the optical system shown in Fig. 1, consisting of one doublet, one aperture, two lenses and an image plane. Let the elements (including the image plane) be labeled from j = 1 to ${\textrm{j}} = 5$ and the ten boundary surfaces be denoted sequentially as i = 1 to ${\textrm{i}} = 10 \equiv {\textrm{n}}$. As shown, the source ray, ${\overline {\textrm{R}} _0} = {\left[ {\begin{array}{cc} {{{\bar{\textrm{P}}}_0}}&{{{\bar{\ell }}_0}} \end{array}} \right]^{\textrm{T}}}$, originates at point source
(1)$${\overline {\textrm{P}} _0} = {\left[ {\begin{array}{cccc} {{\textrm{P}_{0\textrm{x}}}}&{{\textrm{P}_{0\textrm{y}}}}&{{\textrm{P}_{0\textrm{z}}}}&1 \end{array}} \right]^{\textrm{T}}}$$
and travels along the unit directional vector (2)$${\overline {{\ell }}} {}_0 = {\left[ {\begin{array}{cccc} {{\textrm{S}}{\alpha_0}{\textrm{C}}{\beta_0}}&{{\textrm{S}}{\beta_0}}&{{\textrm{C}}{\alpha_0}{\textrm{C}}{\beta_0}}&0 \end{array}} \right]^{\textrm{T}}} \equiv {\left[ {\begin{array}{cccc} {{\ell_{0{\textrm{x}}}}}&{{\ell_{0{\textrm{y}}}}}&{{\ell_{0{\textrm{z}}}}}&0 \end{array}} \right]^{\textrm{T}}}, $$
where S and C denote sine and cosine, respectively. To illustrate the methodology proposed in this study using the system shown in Fig. 1, it is first necessary to determine the position and orientation (i.e., pose) matrices of all the elements (j = 1 to j = 5) and boundary surfaces (i = 1 to i = 10) in the system with respect to the world coordinate frame ${({\textrm{xyz}})_0}$, as described in the following example.Example 1: To perform a tolerance analysis in the design stage of an optical system, each element (say jth element) of the system should possess six degrees-of-freedom (i.e., three position and three rotation variables:${\textrm{t}_{\textrm{ejx}}}$, ${\textrm{t}_{\textrm{ejy}}}$, separation ${\textrm{v}_{\textrm{j}}}$, ${\omega _{\textrm{ejx}}}$, ${\omega _{\textrm{ejy}}}$, ${\omega _{\textrm{ejz}}}$). If the coordinate frame ${({\textrm{xyz}})_{\textrm{ej}}}$ (j = 1 to j = 5) of each element is defined in such a way that its origin is located at the center of its first boundary surface, then the following pose matrices ${}^\textrm{0}{\bar{{\textrm {A}}}_{\textrm{ej}}}$ are obtained for the five elements in Fig. 1:
(3a)$${}^\textrm{0}{\bar{{\textrm {A}}}_{\textrm{e1}}} = \textrm{tran}({{\textrm{t}}_{\textrm{e1x}}}{,}{{\textrm{t}}_{\textrm{e1y}}}{,}{{\textrm{R}}_{1}}){\textrm{rot}}(\bar{{\textrm{z}}},{\omega _{\textrm{e1z}}})\textrm{rot}(\bar{\textrm{y}},{\omega _{\textrm{e1y}}}){\textrm{rot}}(\bar{{\textrm{x}}},{\omega _{\textrm{e1x}}}), $$
(3b)$${}^{0}{\bar{{\textrm {A}}}_{\textrm{e2}}} = {\textrm{tran}}({\textrm{t}}_{\textrm{e2x}}\textrm{,}{\textrm{t}}_{\textrm{e2y}}\textrm{,}{{\textrm{q}}_{\textrm{e1}}} + {{\textrm{q}}^{\prime}_{\textrm{e1}}} + {\textrm{v}_2}){\textrm{rot}(\bar{\textrm {z}},}{\omega _{\textrm{e2z}}})\textrm{rot}(\bar{\textrm{y}},{\omega _{\textrm{e2y}}})\textrm{rot}(\bar{{\textrm{x}}},{\omega _{\textrm{e2x}}}),$$
(3c)$${}^\textrm{0}{\bar{{\textrm {A}}}_{\textrm{e3}}} = \textrm{tran}({\textrm{t}_{\textrm{e3x}}}\textrm{,}{\textrm{t}_{\textrm{e3y}}},{{\textrm{q}}_{\textrm{e1}}} + {{\textrm{q}}^{\prime}_{\textrm{e1}}} + {{\textrm{v}}_2} + {{\textrm{q}}_{\textrm{e2}}} + {{\textrm{v}}_3} + {{\textrm{R}}_6})\textrm{rot}(\bar{\textrm{z}},{\omega _{\textrm{e3z}}})\textrm{rot}(\bar{\textrm{y}},{\omega _{\textrm{e3y}}})\textrm{rot}(\bar{{\textrm{x}}},{\omega _{\textrm{e3x}}}),$$
(3d)$$\begin{aligned}{}^\textrm{0}{{\bar{{\textrm {A}}}}_{\textrm{e4}}} &= \textrm{tran(}{\textrm{t}_{\textrm{e4x}}}\textrm{,}{\textrm{t}_{\textrm{e4y}}}\textrm{,}{{\textrm {q}}_{\textrm{e1}}} + {{{\textrm{q}}^{\prime}}_{\textrm{e1}}} + {{\textrm{v}}_2} + {{\textrm {q}}_{\textrm{e2}}} + {{\textrm{v}}_3} + {{\textrm {q}}_{\textrm{e3}}} + {{\textrm{v}}_4} + {{\textrm{R}}_8})\\ &\quad\textrm{ }\textrm{rot}(\bar{\textrm{z}},{\omega _{\textrm{e4z}}}){\textrm{rot}}(\bar{\textrm{y}},{\omega _{\textrm{e4y}}})\textrm{rot}(\bar{{\textrm{x}}},{\omega _{\textrm{e4x}}}), \end{aligned}$$
(3e)$$\begin{aligned}{}^\textrm{0}{{\bar{{\textrm {A}}}}_{\textrm{e5}}} &= \textrm{tran(}{\textrm{t}_{\textrm{e5x}}}\textrm{,}{\textrm{t}_{\textrm{e5y}}}\textrm{,}{{\textrm {q}}_{\textrm{e1}}} + {{{\textrm{q}}^{\prime}}_{\textrm{e1}}} + {{\textrm{v}}_2} + {{\textrm {q}}_{\textrm{e2}}} + {{\textrm{v}}_3} + {{\textrm {q}}_{\textrm{e3}}} + {{\textrm{v}}_4} + {{\textrm {q}}_{\textrm{e4}}} + {{\textrm{v}}_5})\\ &\quad\textrm{ }\textrm{rot}(\bar{\textrm{z}},{\omega _{\textrm{e5z}}}){\textrm{rot}}(\bar{\textrm{y}},{\omega _{\textrm{e5y}}})\textrm{rot}(\bar{{\textrm{x}}},{\omega _{\textrm{e5x}}}), \end{aligned}$$
where $\textrm{rot}$ and $\textrm{tran}$ are given in Appendix A. ${{\textrm {q}}_{\textrm{ej}}}$ is the thickness of the jth element and ${{\textrm{R}}_{\textrm{i}}}$ is the radius of the ith (${\textrm{i}} \in \{{1,2,3,6,7,8,9} \}$) boundary surface. The pose matrices of the boundary coordinate frames (i.e., ${{\textrm{(xyz)}}_{\textrm{i}}}$, i = 1 to i = 10) with respect to their corresponding element coordinate frames ${\textrm{(xyz)}}_{\textrm{ej}}$ (j = 1 to j = 5) are given by (4)$$\begin{aligned} {}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_{1}} & = {\textrm{I}_{4 \times 4}},\ {}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_{2}} = \textrm{tran(0},0, - {{\textrm{R}}_{1}} + {\textrm {q}}_{\textrm{e1}} + {{\textrm{R}}_{2}}),{}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_3} = {\textrm{tran(0}},0\textrm{,} - {{\textrm{R}}_1} + {{\textrm {q}}_{\textrm{e1}}} + {{{\textrm{q}}^{\prime}}_{\textrm{e1}}} + {{\textrm{R}}_3}),\\ {}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_{2}} &= {\textrm{tran(0}},0\textrm{,} - {{\textrm{R}}_{1}} + {{\textrm {q}}_{\textrm{e1}}} + {{\textrm{R}}_{2}}),\textrm{ }{}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_3} = \textrm{tran(0},0\textrm{,} - {{\textrm{R}}_1} + {{\textrm {q}}_{\textrm{e1}}} + {{{\textrm{q}}^{\prime}}_{\textrm{e1}}} + {{\textrm{R}}_3}),\\ {}^{\textrm{e2}}{{\bar{{\textrm {A}}}}_4} &= {{\bar{\textrm{I}}}_{4 \times 4}},{}^{\textrm{e2}}{{\bar{{\textrm {A}}}}_5} = \textrm{tran(0},0\textrm{,}{{\textrm {q}}_{e3}}),\textrm{ }{}^{\textrm{e3}}{{\bar{{{\textrm {A}}}}}_6} = {{\bar{\textrm{I}}}_{4 \times 4}},\textrm{ }{}^{\textrm{e3}}{{\bar{\textrm{A}}}_7} = \textrm{tran(0},0\textrm{,} - {{\textrm{R}}_6} + {{\textrm {q}}_\textrm{{e3}}} + {{\textrm{R}}_7}),\\ {}^{\textrm{e4}}{{\bar{\textrm{A}}}_{8}} &= {{\bar{\textrm{I}}}_{4 \times 4}},\textrm{ }{}^{\textrm{e4}}{{\bar{\textrm{A}}}_9} = \textrm{tran(0},0\textrm{,} - {{\textrm{R}}_8} + {{\textrm {q}}_\textrm{{e4}}} + {{\textrm{R}}_9}),\textrm{ }{}^{\textrm{e5}}{{\bar{\textrm{A}}}_{\textrm{10}}} = {{\bar{\textrm{I}}}_{4 \times 4}} \end{aligned}, $$
Example 2: If all of the independent variables of the system shown in Fig. 1 are taken as the system variable vector, then ${\bar{{\textrm{X}}}_{\textrm{sys}}}$ (with dimension ${{\textrm {q}}_{\textrm{sys}}} = 51$) has the form (5)$$\begin{aligned} {{\bar{{\textrm{X}}}}_{\textrm{sys}}} = &\left[ {\begin{array}{ccccc} {{{\textrm{P}}_{\textrm{0x}}}}&{{{\textrm{P}}_{\textrm{0y}}}}&{{{\textrm{P}}_{\textrm{0z}}}}&{{\alpha_0}}&{{\beta_0}} \end{array}} \right.\\ & \begin{array}{ccccc} {\textrm{ }{\xi _{\textrm{air}}}}&{{\xi _{\textrm{e1}}}}&{{{\xi ^{\prime}}_{\textrm{e1}}}}&{{\xi _{\textrm{e3}}}}&{\textrm{ }{\xi _{\textrm{e4}}}} \end{array}\\ &\begin{array}{ccccccc} {\textrm{ }{{\textrm{R}}_1}}&{{{\textrm{R}}_2}}&{{{\textrm{R}}_3}}&{{{\textrm{R}}_6}}&{{{\textrm{R}}_7}}&{{{\textrm{R}}_8}}&{{{\textrm{R}}_9}} \end{array}\\ &\begin{array}{ccccc} {\textrm{ }{{\textrm {q}}_{\textrm{e1}}}}&{{{{\textrm{q}}^{\prime}}_{\textrm{e1}}}}&{{{\textrm {q}}_{\textrm{e2}}}}&{{{\textrm {q}}_{\textrm{e3}}}}&{{{\textrm {q}}_{\textrm{e4}}}} \end{array}\\ &\begin{array}{cccccc} {\textrm{ }{\textrm{t}_{\textrm{e1x}}}}&{{\textrm{t}_{\textrm{e1y}}}}&{\textrm{ }}&{{\omega _{\textrm{e1x}}}}&{{\omega _{\textrm{e1y}}}}&{{\omega _{\textrm{e1z}}}} \end{array}\\ &\begin{array}{cccccc} {\textrm{ }{\textrm{t}_{\textrm{e2x}}}}&{{\textrm{t}_{\textrm{e2y}}}}&{{{\textrm{v}}_2}}&{{\omega _{\textrm{e2x}}}}&{{\omega _{\textrm{e2y}}}}&{{\omega _{\textrm{e2z}}}} \end{array}\\ &\begin{array}{cccccc} {\textrm{ }{\textrm{t}_{\textrm{e3x}}}}&{{\textrm{t}_{\textrm{e3y}}}}&{{{\textrm{v}}_3}}&{{\omega _{\textrm{e3x}}}}&{{\omega _{\textrm{e3y}}}}&{{\omega _{\textrm{e3z}}}} \end{array}\\ &\begin{array}{cccccc} {\textrm{ }{\textrm{t}_{\textrm{e4x}}}}&{{\textrm{t}_{\textrm{e4y}}}}&{{{\textrm{v}}_4}}&{{\omega _{\textrm{e4x}}}}&{{\omega _{\textrm{e4y}}}}&{{\omega _{\textrm{e4z}}}} \end{array}\\ &{\left. {\begin{array}{cccccc} {\textrm{t}_{\textrm{e5x}}}&{{\textrm{t}_{\textrm{e5y}}}}&{{{\textrm{v}}_5}}&{\omega_{\textrm{e5x}}}&{{\omega_{\textrm{e5y}}}}&{{\omega_{\textrm{e5z}}}} \end{array}} \right]}^{\textrm{T}}, \end{aligned}$$
where ${\xi _{\textrm{air}}}$ and ${\xi _{\textrm{ej}}}$ (j = 1 to j = 5) are the refractive indices of air and the jth element, respectively. In [15], the variable vector ${\overline {\textrm{X}} _{\textrm{i}}}$ of a boundary surface was defined in terms of its nine independent variables (Fig. 2), i.e., (6)$${\overline {\textrm{X}} _{\textrm{i}}} = {\left[ {\begin{array}{ccccccccc} {{{\textrm{t}}_{\textrm{ix}}}}&{{{\textrm{t}}_{\textrm{iy}}}}&{{{\textrm{t}}_{\textrm{iz}}}}&{{\omega_{\textrm{ix}}}}&{{\omega_{\textrm{iy}}}}&{{\omega_{\textrm{iz}}}}&{{\xi_{{\textrm{i}} - 1}}}&{{\xi_{\textrm{i}}}}&{{{\textrm{R}}_{\textrm{i}}}} \end{array}} \right]^{\textrm{T}}}, $$
where the first six independent variables correspond to the three translational degrees-of-freedom and three rotational degrees-of-freedoms of the boundary surface; and ${\xi _{{\textrm{i}} - 1}}$ and ${\xi _{\textrm{i}}}$ are the refractive indices of the media on either side of the boundary surface, respectively. In the present study, the system variable vector in Eq. (5) is reformulated by replacing all the angular variables with their corresponding cosine and sine trigonometric functions, i.e., (7)$$\begin{aligned}{{\bar{{\textrm{X}}}}_{\textrm{sys}}} = &\left[ {\begin{array}{ccccccc} {{{\textrm{P}}_{\textrm{0x}}}}&{{{\textrm{P}}_{\textrm{0y}}}}&{{{\textrm{P}}_{\textrm{0z}}}}&{{\textrm{C}}{\alpha_0}}&{{\textrm{S}}{\alpha_0}}&{{\textrm{C}}{\beta_0}}&{{\textrm{S}}{\beta_0}} \end{array}} \right.\\ & \begin{array}{ccccc} {\textrm{ }{\xi _{\textrm{air}}}}&{{\xi _{\textrm{e1}}}}&{{{\xi ^{\prime}}_{\textrm{e1}}}}&{{\xi _{\textrm{e3}}}}&{{\xi _{\textrm{e4}}}} \end{array}\\ &\begin{array}{ccccccc} {\textrm{ }{{\textrm{R}}_1}}&{{{\textrm{R}}_2}}&{{{\textrm{R}}_3}}&{{{\textrm{R}}_6}}&{{{\textrm{R}}_7}}&{{{\textrm{R}}_8}}&{{{\textrm{R}}_9}} \end{array}\\ &\begin{array}{ccccc} {\textrm{ }{{\textrm {q}}_{\textrm{e1}}}}&{{{{\textrm{q}}^{\prime}}_{\textrm{e1}}}}&{{{\textrm {q}}_{\textrm{e2}}}}&{{{\textrm {q}}_{\textrm{e3}}}}&{{{\textrm {q}}_{\textrm{e4}}}} \end{array}\\ &\begin{array}{ccccccccc} {\textrm{t}_{\textrm{e1x}}}&{\textrm{t}_{\textrm{e1y}}}&\quad &{\textrm{C}}{\omega _{\textrm{e1x}}}&{{\textrm{S}}{\omega _{\textrm{e1x}}}}&{{\textrm{C}}{\omega _{\textrm{e1y}}}}&{{\textrm{S}}{\omega _{\textrm{e1y}}}}&{{\textrm{C}}{\omega _{\textrm{e1z}}}}&{{\textrm{S}}{\omega _{\textrm{e1z}}}} \end{array}\\ &\begin{array}{ccccccccc} {\textrm{t}_{\textrm{e2x}}}&{\textrm{t}_{\textrm{e2y}}}&{{\textrm{v}}_2}&{\textrm{C}}{\omega _{\textrm{e2x}}}&{{\textrm{S}}{\omega _{\textrm{e2x}}}}&{{\textrm{C}}{\omega _{\textrm{e2y}}}}&{{\textrm{S}}{\omega _{\textrm{e2y}}}}&{{\textrm{C}}{\omega _{\textrm{e2z}}}}&{{\textrm{S}}{\omega _{\textrm{e2z}}}} \end{array}\\ &\begin{array}{ccccccccc} {\textrm{t}_{\textrm{e3x}}}&{\textrm{t}_{\textrm{e3y}}}&{\textrm{v}}_3&{\textrm{C}}{\omega _{\textrm{e3x}}}&{{\textrm{S}}{\omega _{\textrm{e3x}}}}&{{\textrm{C}}{\omega _{\textrm{e3y}}}}&{{\textrm{S}}{\omega _{\textrm{e3y}}}}&{{\textrm{C}}{\omega _{\textrm{e3z}}}}&{{\textrm{S}}{\omega _{\textrm{e3z}}}} \end{array}\\ &\begin{array}{ccccccccc} {\textrm{t}_{\textrm{e4x}}}&{\textrm{t}_{\textrm{e4y}}}&{{\textrm{v}}_4}&{\textrm{C}}{\omega _{\textrm{e4x}}}&{{\textrm{S}}{\omega _{\textrm{e4x}}}}&{{\textrm{C}}{\omega _{\textrm{e4y}}}}&{{\textrm{S}}{\omega _{\textrm{e4y}}}}&{{\textrm{C}}{\omega _{\textrm{e4z}}}}&{{\textrm{S}}{\omega _{\textrm{e4z}}}} \end{array}\\ &{\left. {\begin{array}{ccccccccc} {\textrm{t}_{\textrm{e5x}}}&{\textrm{t}_{\textrm{e5y}}}&{{\textrm{v}}_5}&{\textrm{C}}{\omega_{\textrm{e5x}}}&{{\textrm{S}}{\omega_{\textrm{e5x}}}}&{{\textrm{C}}{\omega_{\textrm{e5y}}}}&{\textrm{S}}{\omega_{\textrm{e5y}}}&{{\textrm{C}}{\omega_{\textrm{e5z}}}}&{{\textrm{S}}{\omega_{\textrm{e5z}}}} \end{array}} \right]^{\textrm{T}}}. \end{aligned}$$
It is noted that the dimension of ${\bar{{\textrm{X}}}_{\textrm{sys}}}$ is increased from ${{\textrm {q}}_{\textrm{sys}}} = 51$ in Eq. (5) to ${{\textrm {q}}_{\textrm{sys}}} = 68$ in Eq. (7). Furthermore, not all of the components in Eq. (7) are independent. Nonetheless, it still possesses the following three advantages:
(1) The numerical values of the Jacobian matrix of a ray with respect to any component of Eq. (
5) can be easily determined from the Jacobian matrix of the components of Eq. (
7). For example,
${{\partial {{\bar{{\textrm{R}}}}_7}} \mathord{\left/ {\vphantom {{\partial {{\bar{R}}_7}} {\partial {\omega_{e2z}}}}} \right.} {\partial {\omega _{\textrm{e2z}}}}}$ can be obtained as
(8)$$\frac{{\partial {{\bar{{\textrm{R}}}}_7}}}{{\partial {\omega _{\textrm{e2z}}}}} = \frac{{\partial {{\bar{{\textrm{R}}}}_7}}}{{\partial ({\textrm{C}}{\omega _{\textrm{e2z}}})}}( - {\textrm{S}}{\omega _{\textrm{e2z}}}) + \frac{{\partial {{\bar{{\textrm{R}}}}_7}}}{{\partial ({\textrm{S}}{\omega _{\textrm{e2z}}})}}({\textrm{C}}{\omega _{\textrm{e2z}}}). $$
However, it is impossible to obtains ${{\partial {{\bar{\textrm{R}}}}_7}} \mathord{\left/ {\vphantom {{\partial {{\bar{{\textrm{R}}}}_7}} {\partial ({\textrm{C}}{\omega_{\textrm{e2z}}})}}} \right.} {\partial ({\textrm{C}}{\omega _{\textrm{e2z}}})}$ and ${{\partial {{\bar{{\textrm{R}}}}_7}} \mathord{\left/ {\vphantom {{\partial {{\bar{{\textrm{R}}}}_7}} {\partial ({\textrm{S}}{\omega_{\textrm{e2z}}})}}} \right.} {\partial ({\textrm{S}}{\omega _{{\textrm{e2z}}})}}}$ from ${{\partial {{\bar{\textrm{R}}}_7}} \mathord{\left/ {\vphantom {{\partial {{\bar{R}}_7}} {\partial {\omega_{\textrm{e2z}}}}}} \right.} {\partial {\omega _{\textrm{e2z}}}}}$.(2) To perform raytracing in a system having n boundary surfaces, the pose matrix
${}^\textrm{0}{\bar{{\textrm {A}}}_{\textrm{i}}}$ (i = 1-n) of each boundary surface must first be known. For example, in running optical design software such as Zemax or OSLO, the required values must first be entered into a spreadsheet. Mathematically,
${}^\textrm{0}{\bar{{\textrm {A}}}_{\textrm{i}}}$ is defined by the combination of the translation and rotation matrices listed in the Appendix
A. In the present study, the components of
${}^\textrm{0}{\bar{{\textrm {A}}}_{\textrm{i}}}$ are polynomial functions expressed in terms of the sine and cosine functions of the angular variables. Therefore, the determination of
${{\textrm{d}({}^\textrm{0}{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{d({}^\textrm{0}{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {d{{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}} \right.} {\textrm{d}{{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}$ (and its higher order derivatives) involves only polynomial differentiation, i.e., trigonometric function calls are not required. As a result, the determination process is easily implemented in computer code. By contrast, when
${\bar{{\textrm{X}}}_{\textrm{sys}}}$ is defined in the form shown in Eq. (
5), the computation of each
${{\textrm{d}({}^\textrm{0}{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{d({}^\textrm{0}{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {d{{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}} \right.} {\textrm{d}{{\bar{{\textrm{X}}}}_{\textrm{sys}}}}}$ not only requires many trigonometric function calls, but also many division, multiplication and addition operations (see pages 233-243 of [
15]). Consequently, the computation complexity is significantly increased.
(3) In optical systems design, optimization methods play a key role in tuning the system variables so as to maximize the system performance. To improve the performance of such optimization techniques, the design variables of the system should be scaled to have the same order. Notably, the sine and cosine functions have ranges of –1 to 1, and hence their orders are the same as those of the refractive indices of optical materials. As a result, the system variable vector defined in Eq. (
7) is beneficial in optimization methods provided that the following constraint equations are employed:
(9a)$${\bar{\Phi }_{1 \textrm{j}}} = \textrm{C}\omega _{\textrm{ejx}}^{\textrm{ }2} + \textrm{S}\omega _{\textrm{ejx}}^{\textrm{ }2} - 1 = 0,\textrm{ }({\textrm{j} = 1,3,4} )$$
(9b)$${\bar{\Phi }_{2 \textrm{j}}} = \textrm{C}\omega _{\textrm{ejy}}^{\textrm{ }2} + \textrm{S}\omega _{\textrm{ejy}}^{\textrm{ }2} - 1 = 0,\textrm{ }({\textrm{j} = 1,3,4} )$$
(9c)$${\bar{\Phi }_{3 \textrm{j}}} = \textrm{C}\omega _{\textrm{ejz}}^{\textrm{ }2} + \textrm{S}\omega _{\textrm{ejz}}^{\textrm{ }2} - 1 = 0,\textrm{ }({\textrm{j} = 1,3,4} )$$
However, optimization methods treat the sine and cosine functions as two different variables, and therefore these constraint equations may not be satisfied. Nevertheless, this drawback can be overcome by assigning higher weighting factors.
Mathematically, a ray ${\overline {\textrm{R}} _{\textrm{i}}}$ at the ith boundary surface is a recursive function ${\bar{\textrm{R}}_{\textrm{i}}} = {\bar{\textrm{R}}_{\textrm{i}}}({{{\bar{\textrm{R}}}_{\textrm{i} - 1}},{{\bar{{\textrm{X}}}}_{\textrm{i}}}} )$ with the given source ray ${\bar{\textrm R}_0}$. Its Jacobian matrix with respect to the system variable vector, ${\bar{\textrm{X}}_{\textrm{sys}}}$, can be determined via the chain rule as:
(10)$$\begin{aligned}\frac{{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}} = \left[ {\begin{array}{cc} {{{\partial {{\bar{\textrm P}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{P}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}}\\ {{{\partial {{\bar{\ell }}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\ell }}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \end{array}} \right] &= \frac{{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\frac{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}} + \ldots + \frac{{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{R}}}_{\textrm{i} - 1}}}}\frac{{\partial {{\bar{\textrm{R}}}_{\textrm{i} - 1}}}}{{\partial {{\bar{\textrm{R}}}_{\textrm{i} - 2}}}}\ldots \frac{{\partial {{\bar{\textrm{R}}}_2}}}{{\partial {{\bar{\textrm{R}}}_1}}}\frac{{\partial {{\bar{\textrm{R}}}_1}}}{{\partial {{\bar{\textrm{X}}}_{1}}}}\frac{{\partial {{\bar{\textrm{X}}}_{1}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}\\ & + \frac{{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{R}}}_{\textrm{i} - 1}}}}\frac{{\partial {{\bar{\textrm{R}}}_{\textrm{i} - 1}}}}{{\partial {{\bar{\textrm{R}}}_{\textrm{i} - 2}}}}\ldots \frac{{\partial {{\bar{\textrm{R}}}_1}}}{{\partial {{\bar{\textrm{R}}}_0}}}\frac{{\partial {{\bar{\textrm{R}}}_0}}}{{\partial {{\bar{\textrm{X}}}_{0}}}}\frac{{\partial {{\bar{\textrm{X}}}_{0}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}. \end{aligned}$$
The term ${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{R}}}_{i - 1}}}}} \right.} {\partial {{\bar{\textrm{R}}}_{i - 1}}}}$ (given in Chapter 7 of [15]) of Eq. (10) is sensitivity matrix of the current ray ${\bar{\textrm{R}}_{\textrm{i}}}$ with the incoming ray ${\bar{\textrm{R}}_{\textrm{i} - 1}}$. In the present study, both the boundary variable vector, ${\bar{\textrm{X}}_{\textrm{i}}}$, and the system variable vector, ${\bar{\textrm{X}}_{\textrm{sys}}}$, are redefined in order to simplify their computation. Therefore, the terms ${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ of Eq. (10) must also be re-derived. To achieve this, it is first appropriate to review the expressions for ray ${\bar{\textrm{R}}_{\textrm{i}}}$ and boundary variable vector ${\bar{\textrm{X}}_{\textrm{i}}}$, respectively, as described in the following section.3. Raytracing at spherical boundary surface
A spherical boundary surface with radius ${{\textrm{R}}_{\textrm{i}}}$ can be expressed in terms of two angle parameters, ${\alpha _{\textrm{i}}}\textrm{ }(0 \le {\alpha _{\textrm{i}}} < 2\pi )$ and ${\beta _{\textrm{i}}}\textrm{ }({{ - \pi } \mathord{\left/ {\vphantom {{ - \pi } 2}} \right.} 2} \le {\beta _{\textrm{i}}} \le {\pi \mathord{\left/ {\vphantom {\pi 2}} \right.} 2})$, by (Fig. 3)
(11)$$^{\textrm i}{\bar{\textrm r}_{\textrm{i}}} = {\left[ {\begin{array}{cccc} {|{{{\textrm{R}}_{\textrm{i}}}} |{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}}}&{|{{{\textrm{R}}_{\textrm{i}}}} |{\textrm{S}}{\beta_{\textrm{i}}}}&{ - |{{{\textrm{R}}_{\textrm{i}}}} |{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}&1 \end{array}} \right]^{\textrm T}}.$$
Furthermore, its two unit normal vectors is given by (12)$$^{\textrm i}{\bar{\textrm n}_{\textrm{i}}} = {\textrm{s}_{\textrm{i}}}{\left[ {\begin{array}{cccc} {{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}}}&{{\textrm{S}}{\beta_{\textrm{i}}}}&{ - {\textrm{C}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}&0 \end{array}} \right]^{\textrm T}}, $$
where ${\textrm{s}_{\textrm{i}}}$ is set to either + 1 or -1 such that the cosine of the incidence angle, ${\theta_{\textrm{i}}}$, is greater than zero (i.e., ${\textrm{C}}{\theta _{\textrm{i}}} > 0$).Equations (11) and (12) are expressed with respect to the boundary coordinate frame, ${(\textrm{xyz})_{\textrm{i}}}$. However, in commercial software packages, the raytracing are generally built relative to the world coordinate frame, ${(\textrm{xyz})_{\textrm{o}}}$. Consequently, the following pose matrix of ${(\textrm{xyz})_{\textrm{i}}}$ with respect to ${(\textrm{xyz})_{\textrm{o}}}$ is required:
(13)$${}^0{\bar{{\textrm {A}}}_{\textrm{i}}} = \left[ {\begin{array}{cccc} {{\textrm{I}_{\textrm {ix}}}}&{{\textrm{J}_{\textrm{ix}}}}&{{\textrm{K}_{\textrm{ix}}}}&{{\textrm{t}_{\textrm{ix}}}}\\ {{\textrm{I}_{\textrm{iy}}}}&{{\textrm{J}_{\textrm{iy}}}}&{{\textrm{K}_{\textrm{iy}}}}&{{\textrm{t}_{\textrm{iy}}}}\\ {{\textrm{I}_{\textrm{iz}}}}&{{\textrm{J}_{\textrm{iz}}}}&{{\textrm{K}_{\textrm{iz}}}}&{{\textrm{t}_{\textrm{iz}}}}\\ 0&0&0&1 \end{array}} \right]. $$
The unit normal vectors ${\bar{\textrm n}_{\textrm{i}}}$ of the boundary surface with respect to frame ${(\textrm {xyz})_{\textrm{o}}}$ can then be obtained via the transformation (14)$${\bar{\textrm n}_{\textrm{i}}} = \left[ {\begin{array}{c} {{\textrm{n}_{\textrm{ix}}}}\\ {{\textrm{n}_{\textrm{iy}}}}\\ {{\textrm{n}_{\textrm{iz}}}}\\ 0 \end{array}} \right] = {}^0{\bar{{\textrm {A}}}_{\textrm{i}}}{}^{\textrm i}{\bar{\textrm n}_{\textrm{i}}} = {\textrm{s}_{\textrm{i}}}\left[ {\begin{array}{c} {{\textrm{I}_{\textrm{ix}}}{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}} + {\textrm{J}_{\textrm{ix}}}{\textrm{S}}{\beta_{\textrm{i}}} - {\textrm{K}_{\textrm{ix}}}{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}\\ {{\textrm{I}_{\textrm{iy}}}{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}} + {\textrm{J}_{\textrm{iy}}}{\textrm{S}}{\beta_{\textrm{i}}} - {\textrm{K}_{\textrm{iy}}}{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}\\ {{\textrm{I}_{\textrm{iz}}}{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}} + {\textrm{J}_{\textrm{iz}}}{\textrm{S}}{\beta_{\textrm{i}}} - {\textrm{K}_{\textrm{iz}}}{\textrm{C}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}\\ 0 \end{array}} \right]. $$
Figure 3 considers a ray ${\bar{\textrm{R}}_{i - 1}} = {\left[ {\begin{array}{cc} {{{\bar{{\textrm {P}}}}_{\textrm{i - 1}}}}&{{{{\bar{\ell }}}_{\textrm{i - 1}}}} \end{array}} \right]^{\textrm{T}}}$ originating at point ${\bar{\textrm{P}}_{\textrm{i - 1}}} = {\left[ {\begin{array}{cccc} {{\textrm{P}_{\textrm{i - 1x}}}}&{{\textrm{P}_{\textrm{i - 1y}}}}&{{\textrm{P}_{\textrm{i - 1z}}}}&1 \end{array}} \right]^{\textrm{T}}}$ and traveling along unit directional vector ${{\bar{\ell }}_{\textrm{i - 1}}} = {\left[ {\begin{array}{cccc} {{\ell_{\textrm{i - 1x}}}}&{{\ell_{\textrm{i - 1y}}}}&{{\ell_{\textrm{i - 1z}}}}&0 \end{array}} \right]^{\textrm{T}}}$ until it undergoes refraction at this ith boundary surface. The incidence point of the ray on the boundary surface is given by (15)$${\bar{\textrm{P}}_{\textrm{i}}} = {\left[ {\begin{array}{cccc} {{\textrm{P}_{\textrm{ix}}}}&{{\textrm{P}_{\textrm{iy}}}}&{{\textrm{P}_{\textrm{iz}}}}&1 \end{array}} \right]^{\textrm T}} = {\left[ {\begin{array}{cccc} {{\textrm{P}_{\textrm{i - 1x}}} + {\ell_{\textrm{i - 1x}}}{\lambda_{\textrm{ }\textrm{i}}}}&{{\textrm{P}_{\textrm{i - 1y}}} + {\ell_{\textrm{i - 1y}}}{\lambda_{\textrm{ }\textrm{i}}}}&{{\textrm{P}_{\textrm{i - 1z}}} + {\ell_{\textrm{i - 1z}}}{\lambda_{\textrm{ }\textrm{i}}}}&1 \end{array}} \right]^{\textrm T}}, $$
where ${\lambda _{\textrm{ }\textrm{i}}}$ is determined by (16)$${\lambda _{\textrm{ }\textrm{i}}} ={-} {{\textrm D}_{\textrm{i}}} \pm \sqrt {{\textrm D_{\textrm{i}}}^2 - {{\textrm E}_{\textrm{i}}}}$$
with (17)$${\textrm{D}_{\textrm{i}}} = {\ell _{\textrm{i} - \textrm{1x}}}({\textrm{P}_{\textrm{i} - \textrm{1x}}} - {\textrm{t}_{\textrm{ix}}}) + {\ell _{\textrm{i} - \textrm{1y}}}({\textrm{P}_{\textrm{i} - \textrm{1y}}} - {\textrm{t}_{\textrm{iy}}}) + {\ell _{\textrm{i} - \textrm{1z}}}({\textrm{P}_{\textrm{i} - \textrm{1z}}} - {\textrm{t}_{\textrm{iz}}}), $$
(18)$${\textrm{E}_{\textrm{i}}} = {({{\textrm{P}}_{\textrm{i} - 1 \textrm{x}}} - {\textrm{t}_{\textrm{ix}}})^2} + {({{\textrm{P}}_{\textrm{i} - 1 \textrm{y}}} - {\textrm{t}_{\textrm{iy}}})^2} + {({{\textrm{P}}_{\textrm{i} - 1 \textrm {z}}} - {\textrm{t}_{\textrm{iz}}})^2} - {\textrm{R}}_{\textrm{i}}^\textrm{2}. $$
Parameters ${\alpha _{\textrm{i}}}$ and ${\beta _{\textrm{i}}}$ in Eq. (11) are given respectively as (19)$${\alpha _{\textrm{i}}} = {\mathop{\textrm {atan}}\nolimits} \,\textrm{2(} - {\tau _{\textrm{i}}}\textrm{ ,}{\sigma _{\textrm{i}}}),$$
(20)$${\beta _{\textrm{i}}} = \textrm{atan}2({\rho _{\textrm{i}}},\sqrt {\sigma _{\textrm{i}}^\textrm{2} + \tau _{\textrm{i}}^2} )\textrm{ },$$
where (21a)$${\sigma _{\textrm{i}}} = {{\textrm I}_{\textrm{ix}}}({{\textrm{P}_{\textrm{i - 1x}}} + {\ell_{\textrm{i - 1x}}}{\lambda_{\textrm{i}}}} )+ {{\textrm I}_{\textrm{iy}}}({{\textrm{P}_{\textrm{i - 1y}}} + {\ell_{\textrm{i - 1y}}}{\lambda_{\textrm{i}}}} )+ {{\textrm I}_{\textrm{iz}}}({{\textrm{P}_{\textrm{i - 1z}}} + {\ell_{\textrm{i - 1z}}}{\lambda_{\textrm{i}}}} )- ({\textrm{I}_{\textrm ix}}{\textrm{t}_{\textrm ix}} + {{\textrm I}_{\textrm{iy}}}{\textrm{t}_{\textrm {iy}}} + {{\textrm I}_{\textrm {iz}}}{\textrm{t}_{\textrm{iz}}}),$$
(21b)$${\rho _{\textrm{i}}} = {{\textrm J}_{\textrm{ix}}}({{\textrm{P}_{\textrm{i - 1x}}} + {\ell_{\textrm{i - 1x}}}{\lambda_{\textrm{i}}}} )+ {{\textrm J}_{\textrm{iy}}}({{\textrm{P}_{\textrm{i - 1y}}} + {\ell_{\textrm{i - 1y}}}{\lambda_{\textrm{i}}}} )+ {{\textrm J}_{\textrm{iz}}}({{\textrm{P}_{\textrm{i - 1z}}} + {\ell_{\textrm{i - 1z}}}{\lambda_{\textrm{i}}}} )- ({\textrm{J}_{\textrm{ix}}}{\textrm{t}_{\textrm{ix}}} + {{\textrm J}_{\textrm{iy}}}{\textrm{t}_{\textrm{iy}}} + {{\textrm J}_{\textrm{iz}}}{\textrm{t}_{\textrm{iz}}}),$$
(21c)$${\tau _{\textrm{i}}} = {{\textrm K}_{\textrm{ix}}}({{\textrm{P}_{\textrm{i - 1x}}} + {\ell_{\textrm{i - 1x}}}{\lambda_{\textrm{i}}}} )+ {{\textrm K}_{\textrm{iy}}}({{\textrm{P}_{\textrm{i - 1y}}} + {\ell_{\textrm{i - 1y}}}{\lambda_{\textrm{i}}}} )+ {{\textrm K}_{\textrm{iz}}}({{\textrm{P}_{\textrm{i - 1z}}} + {\ell_{\textrm{i - 1z}}}{\lambda_{\textrm{i}}}} )- ({\textrm{K}_{\textrm{ix}}}{\textrm{t}_{\textrm{ix}}} + {{\textrm K}_{\textrm{iy}}}{\textrm{t}_{\textrm{iy}}} + {{\textrm K}_{\textrm{iz}}}{\textrm{t}_{\textrm{iz}}}). $$
Function $\textrm{atan2}$ in Eqs. (19) and (20) returns an arctangent in the range of $- \pi$ to $\pi$.The incidence angle ${\theta _{\textrm{i}}}$ of ray ${\bar{\textrm{R}}_{\textrm{i} - 1}} = {\left[ {\begin{array}{cc} {{{\bar{{\textrm {P}}}}_{\textrm{i - 1}}}}&{{{{\bar{\ell }}}_{\textrm{i - 1}}}} \end{array}} \right]^{\textrm T}}$ on the boundary surface has the form
(22)$${\textrm{C}}{\theta _{\textrm{i}}} ={-} {\overline {\ell } _{\textrm{i - 1}}} \bullet {\bar{\textrm n}_{\textrm{i}}} ={-} ({{\ell_{\textrm{i - 1x}}}{\textrm{n}_{\textrm{ix}}} + {\ell_{\textrm{i - 1y}}}{\textrm{n}_{\textrm {iy}}} + {\ell_{\textrm{i - 1z}}}{\textrm{n}_{\textrm {iz}}}} ).$$
In accordance with Snell’s law, the refraction angle ${\underline{\theta } _{\textrm{i}}}$ between two optical media is given as (23)$${\textrm{S}}{\underline{\theta } _{\textrm{i}}} = ({{{{\xi_{\textrm{i - 1}}}} \mathord{\left/ {\vphantom {{{\xi_{\textrm{i - 1}}}} {{\xi_{\textrm{i}}}}}} \right.} {{\xi_{\textrm{i}}}}}} ){\textrm{S}}{\theta _{\textrm{i}}} \equiv {{\textrm N}_{\textrm{i}}}{\textrm{S}}{\theta _{\textrm{i}}}, $$
where ${\xi _{\textrm{i}}}$ is the refractive index of medium i and ${{{\textrm {N}_{\textrm{i}}} \equiv {\xi _{\textrm{i - 1}}}} \mathord{\left/ {\vphantom {{{N_{\textrm{i}}} \equiv {\xi_{\textrm{i - 1}}}} {{\xi_{\textrm{i}}}}}} \right.} {{\xi _{\textrm{i}}}}}$ is the refractive index of medium i-1 relative to that of medium i. Thus, the refracted unit directional vector, ${{\bar{\ell }}_{\textrm{i}}}$, can be expressed as (24)$${{\bar{\ell }}_{\textrm{i}}} = \left[ {\begin{array}{c} {{\ell_{\textrm{ix}}}}\\ {{\ell_{\textrm{iy}}}}\\ {{\ell_{\textrm{iz}}}}\\ 0 \end{array}} \right] = \left[ {\begin{array}{c} { - {\textrm{n}_{\textrm{ix}}}\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}{\textrm{C}}{\theta_{\textrm{i}}})}^2}} + {{\textrm N}_{\textrm{i}}}({\ell_{\textrm{i - 1x}}} + {\textrm{n}_{\textrm{ix}}}{\textrm{C}}{\theta_{\textrm{i}}})}\\ { - {\textrm{n}_{\textrm{iy}}}\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}{\textrm{C}}{\theta_{\textrm{i}}})}^2}} + {{\textrm N}_{\textrm{i}}}({\ell_{\textrm{i - 1y}}} + {\textrm{n}_{\textrm{iy}}}{\textrm{C}}{\theta_{\textrm{i}}})}\\ { - {\textrm{n}_{\textrm{iz}}}\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}{\textrm{C}}{\theta_{\textrm{i}}})}^2}} + {{\textrm N}_{\textrm{i}}}({\ell_{\textrm{i - 1z}}} + {\textrm{n}_{\textrm{iz}}}{\textrm{C}}{\theta_{\textrm{i}}})}\\ 0 \end{array}} \right], $$
where ${\bar{\textrm{n}}_{\textrm{i}}} = {\left[ {\begin{array}{cccc} {{\textrm{n}_{\textrm{ix}}}}&{{\textrm{n}_{\textrm{iy}}}}&{{\textrm{n}_{\textrm{iz}}}}&0 \end{array}} \right]^{\textrm T}}$ and $\textrm{C}{\theta _{\textrm{i}}}$ are given by Eqs. (14) and (22), respectively.In this study, the boundary variable vector, ${\overline {\textrm{X}}}_{\textrm{i}}$, comprises the twelve components of Eq. (13), the refractive indices ${\xi _{i - 1}}$ and ${\xi _{\textrm{i}}}$, and the radius ${{\textrm{R}}_{\textrm{i}}}$, i.e.,
(25)$${\overline {\textrm{X}} _{\textrm{i}}} = {\left[ {\begin{array}{ccccccccccccccc} {{\textrm{I}_{\textrm {ix}}}}&{{\textrm{I}_{\textrm{iy}}}}&{{\textrm{I}_{\textrm{iz}}}}&{{\textrm{J}_{\textrm{ix}}}}&{{\textrm{J}_{\textrm{iy}}}}&{{\textrm{J}_{\textrm{iz}}}}&{{\textrm{K}_{\textrm{ix}}}}&{{\textrm{K}_{\textrm{iy}}}}&{{\textrm{K}_{\textrm{iz}}}}&{{\textrm{t}_{\textrm{ix}}}}&{{\textrm{t}_{\textrm{iy}}}}&{{\textrm{t}_{\textrm{iz}}}}&{{\xi_{\textrm{i} - 1}}}&{{\xi_{\textrm{i}}}}&{{{\textrm{R}}_{\textrm{i}}}} \end{array}} \right]^{\textrm T}}.$$
As will be noted later in the Appendix B, the term ${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ is required to determine the Jacobian matrix of the ray ${\bar{\textrm R}_{\textrm{i}}}$. Notably, ${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ can be easily obtained and readily implemented in computer code if ${\overline {\textrm{X}} _{\textrm{i}}}$ is defined as shown in Eq. (25). Equation (25) is very different from its previous definition given by Eq. (6). By contrast, if ${\overline {\textrm{X}} _{\textrm{i}}}$ is specified as shown in Eq. (6) (i.e., as in [15]), many trigonometric function calls, multiplications and additions are required to determine ${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$.4. Ray Jacobian matrix ${{\partial {{\bar{\textrm R}}_{\textrm{i}}}} / {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$
The ray ${\bar{\textrm{R}}_{\textrm{i}}}$ includes its incidence point ${\bar{\textrm P}_{\textrm i}}$ at the ith boundary surface and its unit directional vector ${\bar{\ell }_{\textrm{i}}}$, i.e., ${\bar{\textrm{R}}_{\textrm{i}}} = {\left[ {\begin{array}{cc} {{{\bar{{\textrm {P}}}}_{\textrm{i}}}}&{{{{\bar{\ell }}}_{\textrm{i}}}} \end{array}} \right]^{\textrm T}}$. Its Jacobian matrix with respect to system variable is determined by Eq. (10), indicating ${{\partial {{\bar{\textrm R}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{R}}_{\textrm{i}}}} {\partial {{\bar{X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$ is required. ${{\partial {{\bar{\textrm R}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{R}}_{\textrm{i}}}} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$ includes two parts, namely the Jacobian matrices ${{\partial {{\bar{\textrm P}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{P}}_{\textrm{i}}}} {\partial {{\bar{X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$ and ${{\partial {{\bar{\ell }}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\ell }}_{\textrm{i}}}} {\partial {{\bar{X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$. The former term can be obtained by differentiating ${\overline {\textrm P} _{\textrm{i}}}$ given in Eq. (15) with respect to the boundary variable vector ${\overline {\textrm{X}} _{\textrm{i}}}$. Since ${{\partial {{\bar{\textrm P}}_{\textrm{i} - 1}}} \mathord{\left/ {\vphantom {{\partial {{\bar{P}}_{i - 1}}} {\partial {{\bar{X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}} = {{\partial {{\bar{\ell }}_{\textrm{i} - 1}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\ell }}_{\textrm{i} - 1}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \bar{0}$, it follows that
(26)$$\frac{{\partial {{\bar{\textrm{P}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = {\left[ {\begin{array}{c} {{{\partial {{\textrm{P}}_{\textrm {ix}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{P}}_{\textrm {ix}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {{\textrm{P}}_{\textrm {iy}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{P}}_{\textrm{iy}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {{\textrm{P}}_{\textrm{iz}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{P}}_{\textrm{iz}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {\bar{0}} \end{array}} \right]_{4 \times 15}} = {\bar{\ell }_{\textrm{i} - 1}}\textrm{ }\frac{{\partial {\lambda _{\textrm{ }\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}. $$
The term ${{\partial {\lambda _{\textrm{ }\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\lambda_{\textrm{ }\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ in Eq. (26) is obtained from Eq. (16) as (27)$$\frac{{\partial {\lambda _{\textrm{ }\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} ={-} \frac{{\partial {{\textrm D}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} \pm \frac{{{\textrm{D}_{\textrm{i}}}}}{{\sqrt {(\textrm{D}_{\textrm{i}}^\textrm{2} - {{\textrm E}_{\textrm{i}}})} }}\frac{{\partial {{\textrm D}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} \pm \frac{{ - 1}}{{2\sqrt {(\textrm{D}_{\textrm{i}}^\textrm{2} - {{\textrm E}_{\textrm{i}}})} }}\frac{{\partial {{\textrm E}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}, $$
in which ${\textrm{D}_{\textrm{i}}}$ and ${\textrm{E}_{\textrm{i}}}$ are given by Eqs. (17) and (18), respectively, and (28)$$\frac{{\partial {\textrm{D}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0&0&0&0&0&{ - {\ell_{\textrm{i - 1x}}}}&{ - {\ell_{\textrm{i - 1y}}}}&{ - {\ell_{\textrm{i - 1z}}}}&0&0&0 \end{array}} \right], $$
(29)$$\frac{{\partial {\textrm{E}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \left[ 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ {2({\textrm{t}_{\textrm{ix}}} - {\textrm{P}_{\textrm{i - 1x}}})}\ \ {2({\textrm{t}_{\textrm{iy}}} - {\textrm{P}_{\textrm{i - 1y}}})}\ \ {2({\textrm{t}_{\textrm{iz}}} - {\textrm{P}_{\textrm{i - 1z}}})}\ \ 0\ \ 0\ \ { - 2{{\textrm{R}}_{\textrm{i}}}} \right].$$
The second part of the ray Jacobian matrix ${{\partial {{\bar{\textrm R}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{R}}_{\textrm{i}}}} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$, i.e., the Jacobian matrix ${{\partial {{\bar{\ell }}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\ell }}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ of the unit directional vector of the refracted ray, can be obtained by differentiating Eq. (24). Since ${{\partial {{\bar{\ell }}_{\textrm{i} - 1}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\ell }}_{\textrm{i} - 1}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}} = \bar{0}$, ${{\partial {{{\bar{\ell }}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{{\bar{\ell }}}_{\textrm{i}}}} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}$ is obtained as (30)$$\begin{aligned}\frac{{\partial {{\bar{\ell }}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = {\left[ {\begin{array}{c} {{{\partial {\ell_{\textrm {ix}}}} \mathord{\left/ {\vphantom {{\partial {\ell_{\textrm{ix}}}} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm X}}_{\textrm{i}}}}}}\\ {{{\partial {\ell_{\textrm {iy}}}} \mathord{\left/ {\vphantom {{\partial {\ell_{\textrm{iy}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {\ell_{\textrm {iz}}}} \mathord{\left/ {\vphantom {{\partial {\ell_{\textrm{iz}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {\bar{0}} \end{array}} \right]_{4 \times 15}} = \left( { - \sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}{\textrm{C}}{\theta_{\textrm{i}}})}^2}} } \right)\textrm{ }\frac{{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\textrm{ + }\left( {\frac{{ - \textrm{N}_{\textrm{i}}^\textrm{2}{\textrm{C}}{\theta_{\textrm{i}}}}}{{\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}{\textrm{C}}{\theta_{\textrm{i}}})}^2}} }}} \right)\textrm{ }\frac{{\partial ({\textrm{C}}{\theta _{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\textrm{ }{{\bar{\textrm{n}}}_{\textrm{i}}}\\ \textrm{ } + \left( {\frac{{{\textrm{N}_{\textrm{i}}}({1 - {\textrm{C}^2}{\theta_{\textrm{i}}}} )}}{{\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}{\textrm{C}}{\theta_{\textrm{i}}})}^2}} }}} \right)\textrm{ }\frac{{\partial {{\textrm N}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\textrm{ }{{\bar{\textrm{n}}}_{\textrm{i}}} + \frac{{\partial {{\textrm N}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\textrm{ }({{{\bar{\ell }}_{\textrm{i - 1}}} + ({\textrm{C}}{\theta_{\textrm{i}}}) {{\bar{\textrm{n}}}_{\textrm{i}}}} )+ {{\textrm N}_{\textrm{i}}}\left( {\frac{{\partial ({\textrm{C}}{\theta_{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\textrm{ }{{\bar{\textrm{n}}}_{\textrm{i}}} + ({\textrm{C}}{\theta_{\textrm{i}}})\frac{{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right). \end{aligned}$$
The terms ${\bar{\textrm{n}}_{\textrm{i}}}$ and ${\textrm{C}}{\theta _{\textrm{i}}}$ in Eq. (30) are given by Eqs. (14) and (22), respectively. The Appendix B describes the method employed to determine term ${{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ in Eq. (30). Furthermore, ${{\partial {{\textrm{N}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{N}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ is determined by differentiating ${{\textrm{N}}_\textrm{{i}}} = {{{\xi _\textrm{i - 1}}} \mathord{\left/ {\vphantom {{{\xi_\textrm{i - 1}}} {{\xi_{\textrm{i}}}}}} \right.} {{\xi _{\textrm{i}}}}}$ with respect to ${\bar{\textrm{X}}_{\textrm{i}}}$, to give (31)$$\frac{{\partial {{\rm{N}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0&0&0&0&0&0&0&0&{{1 \mathord{\left/ {\vphantom {1 {{\xi_{\textrm{i}}}}}} \right.} {{\xi_{\textrm{i}}}}}}&{{{ - {{\rm{N}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{ - {{\rm{N}}_{\textrm{i}}}} {{\xi_{\textrm{i}}}}}} \right.} {{\xi_{\textrm{i}}}}}}&0 \end{array}} \right].$$
Finally, the term ${{\partial ({\textrm{C}}{\theta _{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({\textrm{C}}{\theta_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{{\rm{X}}}}_{\textrm{i}}}}}$ in Eq. (30) can be computed directly from Eq. (22) as (32)$$\frac{{\partial ({\textrm{C}}{\theta _{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} ={-} \left( {{\ell_{\textrm{i - 1x}}}\frac{{\partial {{\rm{n}}_{\textrm{ix}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} + {\ell_{\textrm{i - 1y}}}\frac{{\partial {{\rm{n}}_{\textrm{iy}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} + {\ell_{\textrm{i - 1z}}}\frac{{\partial {{\rm{n}}_{\textrm{iz}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right). $$
Combining Eqs. (26) and (30), the Jacobian matrix ${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ required for Eq. (10) is obtained. Once having ${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ in Eq. (32), we can determine the ray Jacobian matrix with respect to system variable vector by $({{{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} )({{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} )$, where ${{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ converts the boundary variable vector ${\bar{\textrm{X}}_{\textrm{i}}}$ to the system variable vector ${\bar{\textrm{X}}_{\textrm{sys}}}$ (Chapter 8 of [15]). It is notable that ${\bar{\textrm{X}}_{\textrm{i}}}$ is redefined by Eq. (25) in this study, leading the computation of ${{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ only involving polynomial differentiation without trigonometric function calls. Three examples are given in Appendix C to illustrate this fact.This following example will evaluate the computation speed of the methodology proposed in this study for computing the ray Jacobian matrix.
Example 3: Table 1 shows the computation times required by the methods proposed in [15] and in the present study, respectively, to perform the raytracing of ${\bar{\textrm{R}}_{\textrm{i}}}$ and the computation of four Jacobian matrices (i.e.,${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{R}}}_{{\rm{i}} - 1}}}}} \right.} {\partial {{\bar{\textrm{R}}}_{{\rm{i}} - 1}}}}$, ${{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$, ${{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}}}$and ${{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$, i = 1 to i = 9) for the optical system shown in Fig. 1. For comparison purposes, the relative computation speeds are also evaluated by dividing the respective CPU times by that required for tracing ray ${\bar{\textrm{R}}_{\textrm{i}}}$. The computations were performed on a personal computer programmed using FORTRAN with double precision. In performing the computations, the variable values of the optical system were set as follows:
(33)$$\begin{array}{l} {\xi _{{\rm{air}}}} = 1,\textrm{ }{\xi _{{\rm{e}}1}} = 1.65,\textrm{ }{{\xi ^{\prime}}_{{\rm{e}}1}} = 1.71736,\textrm{ }{\xi _{{\rm{e}}3}} = 1.52583,{\xi _{{\rm{e}}4}} = 1.65,\textrm{ }{{\textrm{R}}_1} = 38.2219,\\ {{\textrm{R}}_2} ={-} 56.0857,\textrm{ }{{\textrm{R}}_3} ={-} 590.682,\textrm{ }{{\textrm{R}}_6} ={-} 41.7957,\textrm{ }{{\textrm{R}}_7} = 29.3446,\textrm{ }{{\textrm{R}}_8} = 63.5635,\\ {{\textrm{R}}_9} ={-} 56.8655,\textrm{ }{{\textrm {q}}_{{\rm{e}}1}} = 15.8496,\textrm{ }{{{\textrm{q}}^{\prime}}_{{\rm{e}}1}} = 5.969,\textrm{ }{{\textrm {q}}_{{\rm{e}}3}} = 2.5146,\textrm{ }{{\textrm {q}}_{{\rm{e}}4}} = 6.096,\\ {{\rm{v}}_2} = 3.0226,\textrm{ }{{\rm{v}}_3} = 14.028,{{\rm{v}}_4} = 7.9248,\textrm{ }{{\rm{v}}_5} = 49.6316,\textrm{ }{\omega _{{\rm{e}}1{\rm{x}}}} ={-} {0.2^ \circ },\textrm{ }{\omega _{{\rm{e}}1{\rm{y}}}} ={-} {0.5^ \circ },\\ {\omega _{{\rm{e}}3{\rm{x}}}} = {0.5^ \circ },\textrm{ }{\omega _{{\rm{e}}3{\rm{y}}}} = {1.2^ \circ },\textrm{ }{\omega _{{\rm{e}}4{\rm{x}}}} ={-} {1.2^ \circ },\textrm{ }{\omega _{{\rm{e}}4{\rm{y}}}} ={-} {1^ \circ },\\ {{\textrm {q}}_{{\rm{e2}}}} = {{\rm{t}}_{\rm{{ejx}}}} = {{\rm{t}}_{{\rm{ejy}}}} = {\omega _{{\rm{ejz}}}} = 0,({\rm{j}} = 1 - 5)\\ {\omega _{{\rm{e}}2{\rm{x}}}} = {\omega _{{\rm{e}}2{\rm{y}}}} = {\omega _{{\rm{e}}5{\rm{x}}}} = {\omega _{{\rm{e}}5{\rm{y}}}} = 0 \end{array}$$
As shown in the third row of Table 1, the CPU times required to calculate ${{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ for ${\bar{\textrm{X}}_{\textrm{sys}}}$ defined in Eqs. (5) and (7) are equal to $167.160 \times {10^{ - 5}}$sec and $18.860 \times {10^{ - 5}}$sec, respectively. Moreover, the relative CPU times of the two methods are 66 and 7, respectively. In other words, the computation speed of the proposed method is around ten times faster than that of the method in [15]. The superior performance of the proposed method arises since the computation of ${{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ using Eq. (7) avoids the numerous function calls and other mathematical operations required to compute ${\bar{\textrm{X}}_{\textrm{sys}}}$ using Eq. (5).Table 1. Comparison of CPU times (in ${10^{ - 6}}\textrm{ sec}$) for previous method [15] and current method for tracing of ray $\overline{\textrm{R}}_\textrm{i}$ and computing four Jacobian matrices.
Table 1 also compares the CPU times of the two methods for determining the Jacobian matrixes ${{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}}}$and ${{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ (i = 1 to i = 9), respectively. It is seen that while the CPU time required by the method proposed in this study to compute ${{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}}}$ is longer than that required by the method in [15], the time required to compute ${{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$is reduced. Nonetheless, the performance advantage of the proposed method is degraded by the higher dimension of ${\bar{\textrm{X}}_{\textrm{sys}}}$(i.e., ${{\textrm {q}}_{\textrm{sys}}} = 68$ in the proposed method compared with ${{\textrm {q}}_{\textrm{sys}}} = 51$ in the method proposed in [15]), which increases the upper bound on the number of DO loops performed in the iteration process prior to termination. However, most interested optical systems are axis-symmetrical, leading the dimension ${{\textrm {q}}_{\textrm{sys}}}$ of ${\bar{\textrm{X}}_{\textrm{sys}}}$ is significantly reduce. The advantage of this approach in computation process and coding the computer programs are significant for geometrical optics.
5. Conclusions
The ray Jacobian matrix expresses the first-order derivative of a skew-ray with respect to the system variable vector. However, since the rays in geometrical optics are a highly-nested recursive function of the system variables, computing the Jacobian matrix is extremely complex and has thus been seldom attempted in previous papers. In a previous study by the present group [15], the independent variables of a 3-D optical system were taken as the system variable vector, ${\bar{\textrm{X}}_{\textrm{sys}}}$. However, the method proposed in [15] for determining ${{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{ej}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{ej}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ (and ${{{\textrm{d}}({}^{{\rm{ej}}}{{\bar{\textrm{A}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^{{\rm{ej}}}{{\bar{\textrm{A}}}_{\textrm{i}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$) strictly required the pose matrix ${}^\textrm{0}{\bar{\textrm{A}}_{{\rm{ej}}}}$ (and ${}^{{\rm{ej}}}{\bar{\textrm{A}}_{\textrm{i}}}$) of an element (and a boundary surface) to be specified as a translation motion followed by roll-pitch-yaw motion. To overcome this constraint, the method proposed in the present study for computing the ray Jacobian matrix replaces the angular variables in the original system variable vector, ${\bar{\textrm{X}}_{\textrm{sys}}}$, with their corresponding cosine and sine functions. The variable vector ${\bar{\textrm{X}}_{\textrm{i}}}$ of a boundary surface is redefined in Eq. (25). The illustrative examples presented in this study have shown that the proposed methodology has three important advantages over that reported in [15]: 1) it is valid for any pose matrix specified by an arbitrary combination of rotation and translation motions; 2) the computation process involves only polynomial differentiation (i.e., no trigonometric function calls are required); and 3) the computation speed is significantly improved. As a result, the proposed method can be more readily integrated with existing commercial optical software packages for optical systems design and analysis purposes.
Appendix A:
(34a)$$\textrm{rot}(\bar{\textrm{x}},\theta ) = \left[ {\begin{array}{cccc} 1&0&0&0\\ 0&{\textrm{C}\theta }&{ - \textrm{S}\theta }&0\\ 0&{\textrm{S}\theta }&{\textrm{C}\theta }&0\\ 0&0&0&1 \end{array}} \right]$$
(34b)$$\textrm{rot}(\bar{\textrm{y}},\theta ) = \left[ {\begin{array}{cccc} {\textrm{C}\theta }&0&{\textrm{S}\theta }&0\\ 0&1&0&0\\ { - \textrm{S}\theta }&0&{\textrm{C}\theta }&0\\ 0&0&0&1 \end{array}} \right]$$
(34c)$$\textrm{rot}(\bar{\textrm{z}},\theta ) = \left[ {\begin{array}{cccc} {\textrm{C}\theta }&{ - \textrm{S}\theta }&0&0\\ {\textrm{S}\theta }&{\textrm{C}\theta }&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right]$$
(34d)$$\textrm{tran}({{\textrm{t}}_{\textrm{x}}},{{\textrm{t}}_{\textrm{y}}},{{\textrm{t}}_{\textrm{z}}}) = \left[ {\begin{array}{cccc} 1&0&0&{{{\textrm{t}}_{\textrm{x}}}}\\ 0&1&0&{{{\textrm{t}}_{\textrm{y}}}}\\ 0&0&1&{{{\textrm{t}}_{\textrm{z}}}}\\ 0&0&0&1 \end{array}} \right]$$
(34e)$$\textrm{RPY(}{\theta _{\textrm{z}}},{\theta _{\rm{y}}},{\theta _{\rm{x}}}) = \textrm{rot}(\bar{\textrm{z}},{\theta _{\textrm{z}}})\textrm{rot}(\bar{\textrm{y}},{\theta _{\textrm{y}}})\textrm{rot}(\bar{\textrm{x}},{\theta _{\textrm{x}}})$$
Appendix B:
The difficulty in determining the Jacobian matrix ${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ lies in determining the Jacobian matrix ${{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ of the unit normal vector from Eq. (14). To overcome this difficulty, it is necessary to re-define the boundary variable vector as shown in Eq. (25). Having done so, ${{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ can be obtained by differentiating Eq. (14) with respect to ${\bar{\textrm{X}}_{\textrm{i}}}$, to give
(35)$$\frac{{\partial {{\bar{\textrm{n}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \left[ {\begin{array}{c} {{{\partial {{\textrm{n}}_{\textrm{ix}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{n}}_{\textrm{ix}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {{\textrm{n}}_{\textrm{iy}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{n}}_{\textrm{iy}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {{\textrm{n}}_{\textrm{iz}}}} \mathord{\left/ {\vphantom {{\partial {{\textrm{n}}_{\textrm{iz}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {\bar{0}} \end{array}} \right] = \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\textrm{ }{}^{\rm{i}}{\bar{\textrm{n}}_{\textrm{i}}} + {}^0{\bar{{\textrm {A}}}_{\textrm{i}}}\textrm{ }\frac{{\partial ({}^\textrm{i}{{\bar{\textrm{n}}}_{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}, $$
where
${{\partial ({}^\textrm{i}{{\bar{\textrm{n}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^i{{\bar{\textrm{n}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ is obtained by differentiating Eq. (
12), to give
(36)$$\frac{{\partial ({}^{\rm{i}}{{\bar{\textrm{n}}}_{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = {{\rm{s}}_{\textrm{i}}}\left[ {\begin{array}{c} { - {\textrm{C}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}\\ 0\\ { - {\textrm{C}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}}}\\ 0 \end{array}} \right]\frac{{\partial {\alpha _{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} + {{\rm{s}}_{\textrm{i}}}\left[ {\begin{array}{c} { - {\textrm{S}}{\beta_{\textrm{i}}}{\textrm{C}}{\alpha_{\textrm{i}}}}\\ {{\textrm{C}}{\beta_{\textrm{i}}}}\\ {{\textrm{S}}{\beta_{\textrm{i}}}{\textrm{S}}{\alpha_{\textrm{i}}}}\\ 0 \end{array}} \right]\frac{{\partial {\beta _{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}. $$
It is noted from Eqs. (
35) and (
36) that
${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$,
${{\partial {\alpha _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\alpha_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ and
${{\partial {\beta _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\beta_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ are required, as addressed in the following:
(1) ${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$:
The components of ${{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$, which is a $4 \times 4 \times 15$ matrix, can be determined by differentiating ${}^0{\bar{{\textrm {A}}}_{\textrm{i}}}$ given in Eq. (13) with respect to boundary variable vector ${\overline {\textrm{X}} _{\textrm{i}}}$ given in Eq. (25), to give (see Fig. 4)
(37)$$\begin{array}{l} \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{I}}_{\textrm{ix}}}}} = \left[ {\begin{array}{cccc} 1&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{I}}_{\textrm{iy}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 1&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{I}}_{\textrm{iz}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0\\ 1&0&0&0\\ 0&0&0&0 \end{array}} \right],\\ \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{J}}_{\textrm{ix}}}}} = \left[ {\begin{array}{cccc} 0&1&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{J}}_{\textrm{iy}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&1&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{J}}_{\textrm{iz}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0\\ 0&1&0&0\\ 0&0&0&0 \end{array}} \right],\\ \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{K}}_{\textrm{ix}}}}} = \left[ {\begin{array}{cccc} 0&0&1&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{K}}_{\textrm{iy}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&1&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{K}}_{\textrm{iz}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0\\ 0&0&1&0\\ 0&0&0&0 \end{array}} \right],\\ \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{t}}_{\textrm{ix}}}}} = \left[ {\begin{array}{cccc} 0&0&0&1\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{t}}_{\textrm{iy}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&0&1\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],\frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{t}}_{\textrm{iz}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0\\ 0&0&0&1\\ 0&0&0&0 \end{array}} \right],\\ \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {\xi _{\textrm{i - 1}}}}} = \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {\xi _{\textrm{i}}}}} = \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\textrm{R}}_{\textrm{i}}}}} = \left[ {\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right]. \end{array}$$
(2) ${{\partial {\alpha _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\alpha_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ and ${{\partial {\beta _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\beta_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$:
${{\partial {\alpha _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\alpha_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ and ${{\partial {\beta _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\beta_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ can be obtained by differentiating Eqs. (19) and (20), respectively, with respect to ${\bar{\textrm{X}}_{\textrm{i}}}$, to give
(38)$$\frac{{\partial {\alpha _{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \frac{1}{{\sigma _{\textrm{i}}^\textrm{2} + \tau _{\textrm{i}}^\textrm{2}}}\left( {{\tau_{\textrm{i}}}\frac{{\partial {\sigma_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} - {\sigma_{\textrm{i}}}\frac{{\partial {\tau_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right), $$
(39)$$\frac{{\partial {\beta _{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \frac{{\sqrt {({\sigma_{\textrm{i}}^\textrm{2} + \tau_{\textrm{i}}^2} )} }}{{(\sigma _{\textrm{i}}^\textrm{2} + \rho _{\textrm{i}}^2 + \tau _{\textrm{i}}^2)}}\frac{{\partial {\rho _{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} - \frac{{{\rho _{\textrm{i}}}}}{{(\sigma _{\textrm{i}}^\textrm{2} + \rho _{\textrm{i}}^2 + \tau _{\textrm{i}}^2)\sqrt {({\sigma_{\textrm{i}}^\textrm{2} + \tau_{\textrm{i}}^2} )} }}\left( {{\sigma_{\textrm{i}}}\frac{{\partial {\sigma_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} + {\tau_{\textrm{i}}}\frac{{\partial {\tau_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right). $$
It is noted that
${{\partial {\sigma _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\sigma_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$,
${{\partial {\rho _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\rho_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ and
${{\partial {\tau _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\tau_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ are required in Eqs. (
38) and (
39). All three terms can be obtained by differentiating
(40)$${\bar{\textrm{P}}_{\textrm{i}}} = {}^\textrm{0}{\bar{\textrm{A}}_{\textrm{i}}}\textrm{ }{\left[ {\begin{array}{cccc} {{\sigma_{\textrm{i}}}}&{{\rho_{\textrm{i}}}}&{{\tau_{\textrm{i}}}}&1 \end{array}} \right]^{\textrm{T}}}. $$
with respect to
${\bar{\textrm{X}}_{\textrm{i}}}$, to give
(41)$$\frac{{\partial {{\bar{\textrm{P}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} = \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\left[ {\begin{array}{c} {{\sigma_{\textrm{i}}}}\\ {{\rho_{\textrm{i}}}}\\ {{\tau_{\textrm{i}}}}\\ 1 \end{array}} \right] + {}^0{\bar{{\textrm {A}}}_{\textrm{i}}}\left[ {\begin{array}{c} {{{\partial {\sigma_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\sigma_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {\rho_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\rho_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {\tau_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\tau_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {\bar{0}} \end{array}} \right], $$
where
${{\partial {{\bar{\textrm{P}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{P}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ is given in Eq. (
26). Therefore,
${{\partial {\sigma _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\sigma_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$,
${{\partial {\rho _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\rho_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ and
${{\partial {\tau _{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\tau_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$ can be obtained from Eq. (
41) as
(42)$$\left[ {\begin{array}{c} {{{\partial {\sigma_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\sigma_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}\textrm{ }}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}\textrm{ }}}}\\ {{{\partial {\rho_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\rho_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {{{\partial {\tau_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {\tau_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}}\\ {\bar{0}} \end{array}} \right] = {({{}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}}} )^{ - 1}}\left( {\frac{{\partial {{\bar{\textrm{P}}}_{\textrm{i}}}}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}} - \frac{{\partial ({}^0{{\bar{{\textrm {A}}}}_{\textrm{i}}})}}{{\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}\left[ {\begin{array}{c} {{\sigma_{\textrm{i}}}}\\ {{\rho_{\textrm{i}}}}\\ {{\tau_{\textrm{i}}}}\\ 1 \end{array}} \right]} \right). $$
Appendix C:
To evaluate the ray Jacobian matrix ${{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{\partial {{\bar{\textrm{R}}}_{\textrm{i}}}} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {\partial {{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ from Eq. (10) numerically, one still needs the key Jacobian matrix ${{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{i}}}} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$, which converts the boundary variable vector ${\bar{\textrm{X}}_{\textrm{i}}}$ to the system variable vector ${\bar{\textrm{X}}_{\textrm{sys}}}$. It is noted from the following examples that the computations of $\textrm{ }{{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{\textrm{ej}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{\textrm{ej}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ and ${{{\textrm{d}}({}^{\textrm{ej}}{{\bar{\textrm{A}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^{\textrm{ej}}{{\bar{\textrm{A}}}_{\textrm{i}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}\textrm{ }$are dependent on the manner in which they are defined. It also noted that their computations involve polynomial differentiation only.
Example 4: It is noted from Example 1 that the pose matrix ${}^\textrm{0}{\bar{\textrm{A}}_{\textrm{e1}}}$ of the first element in Fig. 1 is specified using the notation ${\textrm{trans}}({{\textrm{t}}_{\textrm{e1x}}},{{\textrm{t}}_{\textrm{e1y}}},{{\textrm{R}}_1})$ followed by ${\textrm{RPY}}({\omega _{\textrm{e1z}}},{\omega _{\textrm{e1y}}},{\omega _{\textrm{e1x}}})$, i.e.,
(43)$$\begin{array}{l} {}^0{{\bar{\textrm{A}}}_\textrm{e1}} = {\textrm{trans}}({{\textrm{t}}_{\textrm{e1x}}},{{\textrm{t}}_{\textrm{e1y}}},{{\textrm{R}}_1})\textrm{rot}(\bar{\textrm{z},}{\omega _{\textrm{e1z}}})\textrm{rot}(\bar{{\textrm{y}}},{\omega _{\textrm{e1y}}})\textrm{rot}(\bar{\textrm{x},}{\omega _{\textrm{e1x}}})\\ \textrm{ } = \left[ {\begin{array}{cccc} {{\textrm{C}}{\omega_{\textrm{e1y}}}{\textrm{C}}{\omega_{\textrm{e1z}}}}&{{\textrm{S}}{\omega_{\textrm{e1x}}}{\textrm{S}}{\omega_{\textrm{e1y}}}{\textrm{C}}{\omega_{\textrm{e1z}}} - {\textrm{C}}{\omega_{\textrm{e1x}}}{\textrm{S}}{\omega_{\textrm{e1z}}}}&{{\textrm{C}}{\omega_{\textrm{e1x}}}{\textrm{S}}{\omega_{\textrm{e1y}}}{\textrm{C}}{\omega_{\textrm{e1z}}} + {\textrm{S}}{\omega_{\textrm{e1x}}}{\textrm{S}}{\omega_{\textrm{e1z}}}}&{{{\textrm{t}}_{\textrm{e1x}}}}\\ {{\textrm{C}}{\omega_{\textrm{e1y}}}{\textrm{S}}{\omega_{\textrm{e1z}}}}&{{\textrm{C}}{\omega_{\textrm{e1x}}}{\textrm{C}}{\omega_{\textrm{e1z}}} + {\textrm{S}}{\omega_{\textrm{e1x}}}{\textrm{S}}{\omega_{\textrm{e1y}}}{\textrm{S}}{\omega_{\textrm{e1z}}}}&{ - {\textrm{S}}{\omega_{\textrm{e1x}}}{\textrm{C}}{\omega_{\textrm{e1z}}} + {\textrm{C}}{\omega_{{\rm{e1x}}}}{\textrm{S}}{\omega_{{\rm{e1y}}}}{\textrm{S}}{\omega_{{\rm{e1z}}}}}&{{\textrm{t}_{{\rm{e1y}}}}}\\ { - {\textrm{S}}{\omega_{{\rm{e1y}}}}}&{{\textrm{S}}{\omega_{{\rm{e1x}}}}{\textrm{C}}{\omega_{{\rm{e1y}}}}}&{{\textrm{C}}{\omega_{{\rm{e1x}}}}{\textrm{C}}{\omega_{{\rm{e1y}}}}}&{{{\textrm{R}}_1}}\\ 0&0&0&1 \end{array}} \right] \end{array}. $$
When the system variable vector,
${\bar{\textrm{X}}_{\textrm{sys}}}$, is defined in the form shown in Eq. (
7), the following non-zero components of
${{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}} \equiv {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}{({\rm{u}},{\rm{v}},{\rm{w}})_{4 \times 4 \times 68}}$ are obtained:
(44)$$\begin{array}{l} {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,1,29) = {\textrm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,1,31) = {\textrm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,29) = {\textrm{S}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,32) = {\textrm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,1,30) ={-} 1,{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,28) = {\textrm{S}}{\omega _{{\rm{e1y}}}}{\textrm{C}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,30) = {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,31) = {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,27) ={-} {\textrm{S}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,32) ={-} {\textrm{C}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,27) = {\textrm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,31) = {\textrm{C}}{\omega _{{\rm{e1x}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,28) = {\textrm{S}}{\omega _{{\rm{e1y}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,30) = {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,32) = {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,28) = {\textrm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,29) = {\textrm{S}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,3,27) = {\textrm{S}}{\omega _{{\rm{e1y}}}}{\textrm{C}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,3,30) = {\textrm{C}}{\omega _{{\rm{e1x}}}}{\textrm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,3,31) = {\textrm{C}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,3,28) = {\textrm{S}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,3,32) = {\textrm{S}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,28) ={-} {\textrm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,31) ={-} {\textrm{S}}{\omega _{{\rm{e1x}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,27) = {\textrm{S}}{\omega _{{\rm{e1y}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,30) = {\textrm{C}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,32) = {\textrm{C}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},\\ {\textrm{d}}{A_{{\rm{e1}}}}(3,3,27) = {\textrm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,3,29) = {\textrm{C}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,1,25) = 1,{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,2,26) = 1,\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,3,13) = 1. \end{array}$$
In other words, the determination of
$\textrm{ }{{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{ej}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{ej}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ involves only polynomial differentiation without any trigonometric function calls. However, if the system variable vector is defined in the original form shown in Eq. (
5), determining
$\textrm{ }{{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{ej}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{ej}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ poses a significant challenge.
Example 5: Matrix multiplication does not have a commutative property. Hence, the order in which the matrices are specified in defining ${}^\textrm{0}{\bar{\textrm{A}}_{{\rm{ej}}}}$ is extremely important. In practice, the pose matrix ${}^\textrm{0}{\bar{\textrm{A}}_{{\rm{e1}}}}$ can be obtained via various combinations of translation and rotation motions, e.g.,
(45)$$ =\left[\begin{array}{cccc}{\mathrm{C\omega}_{\mathrm{ely}} \mathrm{C\omega}_{\mathrm{elz}}} & {-\mathrm{C\omega}_{\mathrm{ely}} \mathrm{S\omega}_{\mathrm{elz}}} & {\mathrm{S\omega}_{\mathrm{ely}}} & {\mathrm{t}_{\mathrm{elx}}} \\ {\mathrm{S\omega}_{\mathrm{elx}} \mathrm{S} \omega_{\mathrm{ely}} \mathrm{Co}_{\mathrm{elz}}+\mathrm{C} \omega_{\mathrm{elx}} \mathrm{S} \omega_{\mathrm{elz}}} & {\mathrm{C\omega}_{\mathrm{elx}} \mathrm{Co}_{\mathrm{elz}}-\mathrm{S\omega}_{\mathrm{elx}} \mathrm{S} \omega_{\mathrm{elz}}} & {-\mathrm{S\omega}_{\mathrm{elx}} \mathrm{Co}_{\mathrm{ely}}} & {\mathrm{t}_{\mathrm{ely}} \mathrm{C} \omega_{\mathrm{elx}}} \\ {-\mathrm{C} \omega_{\mathrm{elx}} \mathrm{S} \omega_{\mathrm{ely}} \mathrm{C} \omega_{\mathrm{elz}}+\mathrm{S} \omega_{\mathrm{elx}} \mathrm{S} \omega_{\mathrm{elz}}} & {\mathrm{S} \omega_{\mathrm{elx}} \mathrm{C} \omega_{\mathrm{elz}}+\mathrm{C} \omega_{\mathrm{elx}} \mathrm{S} \omega_{\mathrm{ely}} \mathrm{S} \omega_{\mathrm{elz}}} & {\mathrm{C\omega}_{\mathrm{elx}} \mathrm{C} \omega_{\mathrm{ely}}} & {\mathrm{t}_{\mathrm{ely}} \mathrm{S} \omega_{\mathrm{elx}}+\mathrm{R}_{1}} \\ {0} & {0} & {0} & {1}\end{array}\right] $$
The numerical values of
${{\rm{t}}_{{\rm{e1x}}}}$,
${{\rm{t}}_{{\rm{e1y}}}}$,
${\omega _{{\rm{e1x}}}}$,
${\omega _{{\rm{e1y}}}}$ and
${\omega _{{\rm{e1z}}}}$ in Eq. (
45) are, of course, different from those in Eq. (
44). If the system variable vector,
${\bar{\textrm{X}}_{\textrm{sys}}}$, is defined as shown in Eq. (
7), then the following non-zero components of
$\textrm{ }{{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}} \equiv {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}{({\rm{u}},{\rm{v}},{\rm{w}})_{4 \times 4 \times 68}}$ for
${}^\textrm{0}{\bar{\textrm{A}}_{{\rm{e1}}}}$ in Eq. (
45) are obtained:
(46)$$\begin{array}{l} {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,1,29) = {\rm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,1,31) = {\rm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,28) = {\textrm{S}}{\omega _{{\rm{e1y}}}}{\rm{C}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,30) = {\textrm{S}}{\omega _{{\rm{e1x}}}}{\rm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,31) = {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,27) = {\textrm{S}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,1,32) = {\rm{C}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,1,27) ={-} {\textrm{S}}{\omega _{{\rm{e1y}}}}{\rm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,1,30) ={-} {\rm{C}}{\omega _{{\rm{e1x}}}}{\rm{C}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,1,31) ={-} {\rm{C}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,1,28) = {\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,1,32) = {\textrm{S}}{\omega _{{\rm{e1x}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,29) ={-} {\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,2,32) ={-} {\rm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,27) = {\rm{C}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,31) = {\rm{C}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,28) ={-} {\textrm{S}}{\omega _{{\rm{e1y}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,30) ={-} {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,2,32) ={-} {\textrm{S}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,28) = {\rm{C}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,31) = {\textrm{S}}{\omega _{{\rm{e1x}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,27) = {\textrm{S}}{\omega _{{\rm{e1y}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,30) = {\rm{C}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1z}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,2,32) = {\rm{C}}{\omega _{{\rm{e1x}}}}{\textrm{S}}{\omega _{{\rm{e1y}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(1,3,30) = 1,{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,28) ={-} {\rm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(2,3,29) ={-} {\textrm{S}}{\omega _{{\rm{e1x}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,3,27) = {\rm{C}}{\omega _{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(3,3,29) = {\rm{C}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,1,25) = 1,\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,2,26) = {\rm{C}}{\omega _{{\rm{e1x}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,2,27) = {{\rm{t}}_{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,3,26) = {\textrm{S}}{\omega _{{\rm{e1x}}}},\\ {\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,3,28) = {{\rm{t}}_{{\rm{e1y}}}},{\textrm{d}}{{\rm{A}}_{{\rm{e1}}}}(4,3,13) = 1. \end{array}$$
It is noted from Examples 4 and 5 that the expressions for the Jacobian matrix
${{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ depend on the order of the translation and rotation motions. However, the methodology proposed in [
15] for solving
${{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^\textrm{0}{{\bar{\textrm{A}}}_{{\rm{e1}}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ is valid only for the particular case of
${\rm{trans}}({{\rm{t}}_{{\rm{e1x}}}},{{\rm{t}}_{{\rm{e1y}}}},{{\textrm{R}}_1}){\rm{RPY}}({\omega _{{\rm{e1z}}}},{\omega _{{\rm{e1y}}}},{\omega _{{\rm{e1x}}}})$.
Example 6: The determination of$\textrm{ }{{{\textrm{d}}({}^{{\rm{ej}}}{{\bar{\textrm{A}}}_{\textrm{i}}})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^{{\rm{ej}}}{{\bar{\textrm{A}}}_{\textrm{i}}})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}$ for Eq. (4) is comparatively easier since ${}^{{\rm{ej}}}{\bar{\textrm{A}}_{\textrm{i}}}$ is usually a simple matrix. Consider the matrix
(47)$${}^{\textrm{e1}}{\bar{{\textrm {A}}}_3} = \textrm{tran(0},0\textrm{,} - {{\textrm{R}}_1} + {{\textrm {q}}_{{\rm{e1}}}} + {{\textrm{q}}^{\prime}_{{\rm{e1}}}} + {{\textrm{R}}_3}) = \left[ {\begin{array}{cccc} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&{ - {{\textrm{R}}_1} + {{\textrm {q}}_{{\rm{e1}}}} + {{{\textrm{q}}^{\prime}}_{{\rm{e1}}}} + {{\textrm{R}}_3}}\\ 0&0&0&0 \end{array}} \right]$$
for illustration purposes. When the system variable vector,
${\bar{\textrm{X}}_{\textrm{sys}}}$, is defined as shown in Eq. (
7), the following non-zero components of
${{{\textrm{d}}({}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_3})} \mathord{\left/ {\vphantom {{{\textrm{d}}({}^{\textrm{e1}}{{\bar{{\textrm {A}}}}_3})} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}}} \right.} {{\textrm{d}}{{\bar{\textrm{X}}}_{\textrm{sys}}}}} \equiv {\textrm{d}}{}^{\textrm{e1}}{{\rm{A}}_3}{({\rm{u}},{\rm{v}},{\rm{w}})_{4 \times 4 \times 68}}$ are obtained:
${\textrm{d}}{}^{{\rm{e1}}}{{\rm{A}}_3}(3,4,13) ={-} 1$, ${\textrm{d}}{}^{{\rm{e1}}}{{\rm{A}}_3}(3,4,20) = {\textrm{d}}{}^{{\rm{e1}}}{{\rm{A}}_3}(3,4,21) = {\textrm{d}}{}^{{\rm{e1}}}{{\rm{A}}_3}(3,4,15) = 1$.
Funding
Ministry of Science and Technology, Taiwan (MOST 106-2221-E-006-091-MY3).
Disclosures
The author declare no conflicts of interest.
References
1. R. J. Leveque, Finite Difference Methods for Ordinary and Partial Differential Equations, Steady-State and Time-Dependent Problems, (SIAM, 2007), pp. 3–4.
2. B. D. Stone, “Determination of initial ray configurations for asymmetric systems,” J. Opt. Soc. Am. A 14(12), 3415–3429 (1997). [CrossRef]
3. R. Shi and J. Kross, “Differential ray tracing for optical design,” Proc. SPIE 3737, 149–160 (1999). [CrossRef]
4. T. B. Andersen, “Optical aberration functions: derivatives with respect to axial distances for symmetrical systems,” Appl. Opt. 21(10), 1817–1823 (1982). [CrossRef]
5. T. B. Andersen, “Optical aberration functions: derivatives with respect to surface parameters for symmetrical systems,” Appl. Opt. 24(8), 1122–1129 (1985). [CrossRef]
6. D. P. Feder, “Calculation of an optical merit function and its derivatives with respect to the system parameters,” J. Opt. Soc. Am. 47(10), 913–925 (1957). [CrossRef]
7. D. P. Feder, “Differentiation of ray-tracing equations with respect to constructional parameters of rotationally symmetric systems,” J. Opt. Soc. Am. 58(11), 1494–1505 (1968). [CrossRef]
8. Stavroudis, “A simpler derivation of the formulas for generalized ray tracing,” J. Opt. Soc. Am. 66(12), 1330–1333 (1976). [CrossRef]
9. J. Kross, “Differential ray tracing formulae for optical calculations: Principles and applications,” Proc. SPIE 1013, 10–18 (1989). [CrossRef]
10. F. W. Oertmann, “Differential ray tracing formulae; applications especially to aspheric optical systems,” Proc. SPIE 1013, 20–26 (1989). [CrossRef]
11. W. Mandler, “Uber die Berechnung einfacher GauB-Objektive,” Optik 55, 219–240 (1980).
12. B. D. Stone and G. W. Forbes, “Differential ray tracing in inhomogeneous media,” J. Opt. Soc. Am. A 14(10), 2824–2836 (1997). [CrossRef]
13. ZEMAX is a registered trademark of ZEMAX Development Corporation, Bellevue, Washington, USA. http://www.zemax.com/
14. Code-V is a registered trademark of Optical research Associates, Pasadena, California, USA. http://www.opticalres.com/
15. P. D. Lin, Advanced Geometrical Optics (Springer, 2017).
16. P. D. Lin and R. B. Johnson, “Seidel aberration coefficients: an alternative computational method,” Opt. Express 27(14), 19712–19725 (2019). [CrossRef]