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

 

Fig. 1. Illustrative optical system consisting of one doublet, one aperture, two lenses and an image plane.

Download Full Size | PDF

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:

 

Fig. 2. Boundary variables of a boundary surface.

Download Full Size | PDF

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:

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:

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)

 

Fig. 3. Raytracing at spherical boundary surface.

Download Full Size | PDF

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:

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

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.,

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

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:

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.

View Table

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}$$

 

Fig. 4. Schematic representation of matrix ${{\partial \textrm{(}{}^0{{\bar{\textrm{A}}}_{\textrm{i}}})} / {\partial {{\bar{\textrm{X}}}_{\textrm{i}}}}}$with dimension $4 \times 4 \times 15$.

Download Full Size | PDF

(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]