1. Introduction
Geometrical optics is a fundamental science that has attracted the attention of many researchers over the years. Broadly speaking, existing literature on geometrical optics falls into two main categories, namely system design and system analysis. System design involves determining the specific system variable vector required to create an optical system that meets a particular set of needs. By contrast, system analysis involves investigating the performance of a given system while changing the source variable vector and keeping all of the other system parameters unchanged. This study focuses to system analysis to investigate some fundamental problems of geometrical optics based on the Taylor series expansion of a skew ray.
Skew ray tracing equations are recursive functions. Consequently, the skew ray on the image plane of a system with n boundary surfaces (denoted as $\bar{\textrm{R}}_{\textrm{n}}$) is a highly composite function. As a result some previous studies using commercial ray tracing software to investigate problems such as evaluating the point spread function and modulation transfer function (p. 372 of [1]). However, various studies have shown that many problems of geometrical optics can be investigated more easily using ray derivatives (e.g., [2–6]). For example, Burkhard and Shealy proposed a method for determining the principal radii of curvatures on the caustic surfaces [6]. In the present study, the Taylor series expansion is used to expand $\bar{\textrm{R}}_{\textrm{n}}$ with respect to the source ray variable vector ${\bar{\textrm{X}}_0}$ up to the fifth-order. It is shown that the previously derived $6 \times 6$ paraxial ray tracing equations for non-axially symmetrical systems [7] stem from the first-order expansion, ${{{\mathrm{\partial}} {\bar{\textrm{R}}_{\textrm{n}}}}\mathord{\left/ {\vphantom {{{\mathrm{\partial}} {{\bar{\textrm{R}}_{\textrm{n}}}}} {{\mathrm{\partial}} {{\bar{X}}_0}}}} \right. } {{\mathrm{\partial}} {{\bar{\textrm{X}}}_0}}} $. The conventional $2 \times 2$ translation and refraction/reflection matrices for axis-symmetrical systems can thus be obtained from the corresponding $6 \times 6$ translation and refraction/reflection matrices of a non-axially symmetrical system. Moreover, the numerical results obtained for the point spread function, caustic surfaces and modulation transfer function based on ${{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_0}}}} \right.} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_0}}}$ (chap. 13 of [3]) are more accurate than those derived using the ray tracing method (pp. 361-372 of [1]). In addition, the image orientation analysis of prism systems can be performed using the sub-matrix of ${{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_0}}}} \right.} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_0}}}$ without the need for the trial-and-error method described on page 100 of [1].
The earliest report on ray aberrations was published by Philip Ludwig von Seidel in 1856. Since then, many approaches for deriving equations to compute the Seidel primary ray aberration coefficients have been proposed [8–12]. One of the most widely used methods is that developed by Buchdahl [9], which can compute primary ray aberrations as well as in principle any higher-order term. In his approach, the marginal and chief paraxial rays are traced and the results are then used to compute the unconverted third-order Buchdahl aberration coefficients which correspond to the Seidel aberration coefficients in general meaning. These coefficients can also easily be converted into transverse, longitudinal, and wave aberration values (Sec. 4.4 of [8]). Many of the methods proposed in [8–12] for determining the aberration coefficients require significant algebraic effort. However, the present group recently showed that the Seidel primary ray aberration coefficients can be accurately determined by the third-order expansion of ray ${\bar{\textrm{R}}_{\textrm{n}}}$ [13]. This finding suggests that it is possible to determine the higher-order ray aberrations using higher-order expansions of ${\bar{\textrm{R}}_{\textrm{n}}}$ with respect to ${\bar{\textrm{X}}_0}$. Accordingly, the present study uses the Taylor series expansion of the skew ray on the image plane to show the roadmap for various problems of geometrical optics, including the paraxial ray tracing equations, the point spread function, the caustic surfaces, the modulation transfer function, and the secondary and higher order ray aberrations.
2. Ray tracing equations
Consider an optical system consisting of n boundary surfaces. Let the boundary surfaces be denoted sequentially as i=1 to ${\rm i} = {\rm n}$. By convention, label i=0 is assigned to the source ray ${\overline{\textrm{R}}_0} = {\left[ {\begin{array}{cc} {{{\bar{\textrm{P}}}_0}} &{{{\bar{\ell}}_0}} \end{array}} \right]^{\textrm{T}}}$ originating at point source
(1)$${\bar{\textrm{P}}_0} = {\left[ {\begin{array}{ccc} {{\textrm{P}_{0\textrm{x}}}}&{{\textrm{P}_{0\textrm{y}}}}&{{\textrm{P}_{0\textrm{z}}}} \end{array}} \right]^{\textrm{T}}}$$
and travels along the unit directional vector (2)$$\overline{\ell } {}_0 = {\left[ {\begin{array}{ccc} {\textrm{S}{\mathrm{\alpha}_0}\textrm{C}{\mathrm{\beta}_0}}&{\textrm{S}{\mathrm{\beta}_0}}&{\textrm{C}{\mathrm{\alpha}_0}\textrm{C}{\mathrm{\beta}_0}} \end{array}} \right]^{\textrm{T}}} \equiv {\left[ {\begin{array}{ccc} {{\ell_{\textrm{0x}}}}&{{\ell_{\textrm{0y}}}}&{{\ell_{\textrm{0z}}}} \end{array}} \right]^{\textrm{T}}},$$
where S and C denote sine and cosine, respectively. ${{\alpha} _0}$ is the angle between the meridional plane and $\overline{\ell } {}_0$. ${\mathrm{\beta} _0}$ is the angle between $\overline{\ell } {}_0$ and the projection of $\overline{\ell } {}_0$ on the ${\rm x}{}_0{\rm z}{}_0$ plane. If we define (3)$${\bar{\textrm{X}}_{0}} = {\left[ {\begin{array}{ccccc} {{\textrm{P}_{\textrm{0x}}}}&{{\textrm{P}_{\textrm{0y}}}}&{{\textrm{P}_{\textrm{0z}}}}&{{\mathrm{\alpha}_0}}&{{\mathrm{\beta}_0}} \end{array}} \right]^{\textrm{T}}} \equiv [{{\textrm{x}_0}} ]$$
as the vector of independent variables of the source ray, then we have: (4)$${\bar{\textrm{R}}_0} = {\bar{\textrm{R}}_0}({\bar{\textrm{X}}_0}).$$
Equation (4) indicates the source ray ${\bar{\textrm{R}}_0}$ is function of source ray variable vector ${\bar{\textrm{X}}_{0}}$.Figure 1 shows a 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}}}$ originating at point ${\bar{\textrm{P}}_{\textrm{i-1}}} = {\left[ {\begin{array}{ccc} {{\textrm{P}_{\textrm{i-1x}}}}&{{\textrm{P}_{\textrm{i-1y}}}}&{{\textrm{P}_{\textrm{i-1z}}}} \end{array}} \right]^{\textrm{T}}}$ and traveling along the unit directional vector ${{\bar{\ell }}_{\textrm{i-1}}} = {\left[ {\begin{array}{ccc} {{\ell_{\textrm{i-1x}}}}&{{\ell_{\textrm{i-1y}}}}&{{\ell_{\textrm{i-1z}}}} \end{array}} \right]^{\textrm{T}}}$ in 3-D space until it is reflected/refracted at the ith boundary surface, ${\bar{\textrm{r}}_{\textrm{i}}}$. The incidence point ${\bar{\textrm{P}}_{\textrm{i}}}$ on the boundary is given by (pp. 34-39 of [3])
(5)$${\bar{\textrm{P}}_{\textrm{i}}} = \left[ {\begin{array}{c} {{\textrm{P}_{\textrm{ix}}}}\\ {{\textrm{P}_{\textrm{iy}}}}\\ {{\textrm{P}_{\textrm{iz}}}} \end{array}} \right] = \left[ {\begin{array}{c} {{\textrm{P}_{\textrm{i-1x}}} + {\ell_{\textrm{i-1x}}}{\lambda_{\textrm{i}}}}\\ {{\textrm{P}_{\textrm{i-1y}}} + {\ell_{\textrm{i-1y}}}{\lambda_{\textrm{i}}}}\\ {{\textrm{P}_{\textrm{i-1z}}} + {\ell_{\textrm{i-1z}}}{\lambda_{\textrm{i}}}} \end{array}} \right].$$
After noting ${\bar{\textrm{P}}_{\textrm{i}}}$ is a point of the sphere ${\bar{\textrm{r}}_{\textrm{ i}}}$ and ${{\bar{\ell }}_{\textrm{i-1}}}$ is a unit directional vector having $\ell _{\textrm{i-1x}}^{2} + \ell _{\textrm{i-1y}}^{2} + \ell _{\textrm{i-1z}}^{2} = 1$, one can have ${\lambda _{\textrm{i}}}$ from Eq. (5), given as:
(6)$${\lambda _{\textrm{i}}} ={-} {\textrm{D}_{\textrm{i}}} \pm \sqrt {{\textrm{D}_{\textrm{i}}}^2 - {\textrm{E}_{\textrm{i}}}} ,$$
where (7a)$${\textrm{D}_{\textrm{i}}} ={-} {\textrm{t}_{\textrm{ix}}}{\ell _{\textrm{i-1x}}} - {\textrm{t}_{\textrm{iy}}}{\ell _{\textrm{i-1y}}} - {\textrm{t}_{\textrm{iz}}}{\ell _{\textrm{i-1z}}} + {\textrm{P}_{\textrm{i-1x}}}{\ell _{\textrm{i-1x}}} + {\textrm{P}_{\textrm{i-1y}}}{\ell _{\textrm{i-1y}}} + {\textrm{P}_{\textrm{i-1z}}}{\ell _{\textrm{i-1z}}},$$
(7b)$${\textrm{E}_{\textrm{i}}} = \textrm{P}_{\textrm{i-1x}}^{2} + \textrm{P}_{\textrm{i-1y}}^{2} + \textrm{P}_{\textrm{i-1z}}^{2} - \textrm{R}_{\textrm{i}}^{2} + \textrm{t}_{\textrm{ix}}^{2} + \textrm{t}_{\textrm{iy}}^{2} + \textrm{t}_{\textrm{iz}}^{2} - 2({\textrm{t}_{\textrm{ix}}}{\textrm{P}_{\textrm{i-1x}}} + {\textrm{t}_{\textrm{iy}}}{\textrm{P}_{\textrm{i-1y}}} + {\textrm{t}_{\textrm{iz}}}{\textrm{P}_{\textrm{i-1z}}}).$$
Note that ${\textrm{R}_{\textrm{i}}}$ in Eq. (7b) is the radius of the sphere centered at ${\bar{\textrm{o}}_{\textrm{i}}}\textrm{ = (}{\textrm{t}_{\textrm{ix}}},{\textrm{t}_{\textrm{iy}}},{\textrm{t}_{\textrm{iz}}})$ in coordinate frame ${(\textrm{xyz})_0}$. The parameter ${\lambda_{\textrm{i}}}$ represents the geometrical path length from points ${\overline{P}_{\textrm{i-1}}}$ to ${\overline{P}_{\textrm{i}}}$. Also note that the ± sign in Eq. (6) indicates the two possible intersection points of the ray with a complete sphere. The product of ${\lambda _{\textrm{i}}}$ with the refractive index ${{\mathrm{\xi}} _{\textrm{i-1}}}$ (i.e., ${{\mathrm{\xi}} _{\textrm{i-1}}}{\lambda_{\textrm{i}}}$) is the optical path length between points ${\bar{\textrm{P}}_{\textrm{i-1}}}$ and ${\bar{\textrm{P}}_{\textrm{i}}}$. The reflected and refracted unit directional vectors are given respectively as (8a)$${{\bar{\ell }}_{\textrm{i}}} = \left[ {\begin{array}{c} {{\ell_{\textrm{i-1x}}} + 2{\textrm{n}_{\textrm{ix}}}\textrm{C}{\mathrm{\theta}_{\textbf{i}}}}\\ {{\ell_{\textrm{i-1y}}} + 2{\textrm{n}_{\textrm{iy}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}}}\\ {{\ell_{\textrm{i-1z}}} + 2{\textrm{n}_{\textrm{iz}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}}} \end{array}} \right] \equiv \left[ {\begin{array}{c} {{\ell_{\textrm{ix}}}}\\ {{\ell_{\textrm{iy}}}}\\ {{\ell_{\textrm{iz}}}} \end{array}} \right],$$
(8b)$${{\bar{\ell }}_{\textrm{i}}} = \left[ {\begin{array}{c} { - {\textrm{n}_{\textrm{ix}}}\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}})}^2}} + {\textrm{N}_{\textrm{i}}}({\ell_{\textrm{i-1x}}} + {\textrm{n}_{\textrm{ix}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}})}\\ { - {\textrm{n}_{\textrm{iy}}}\sqrt {1 - \textrm{N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}})}^2}} + {\textrm{N}_{\textrm{i}}}({\ell_{\textrm{i-1y}}} + {\textrm{n}_{\textrm{iy}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}})}\\ { - {\textrm{n}_{\textrm{iz}}}\sqrt {\textrm{1 - N}_{\textrm{i}}^\textrm{2} + {{({\textrm{N}_{\textrm{i}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}})}^2}} + {\textrm{N}_{\textrm{i}}}({\ell_{\textrm{i-1z}}} + {\textrm{n}_{\textrm{iz}}}\textrm{C}{\mathrm{\theta}_{\textrm{i}}})} \end{array}} \right] \equiv \left[ {\begin{array}{c} {{\ell_{\textrm{ix}}}}\\ {{\ell_{\textrm{iy}}}}\\ {{\ell_{\textrm{iz}}}} \end{array}} \right],$$
where ${\bar{\textrm{n}}_{\textrm{i}}} = {\left[ {\begin{array}{ccc} {{\textrm{n}_{\textrm{ix}}}}&{{\textrm{n}_{\textrm{iy}}}}&{{\textrm{n}_{\textrm{iz}}}} \end{array}} \right]^{\textrm{T}}}$ is the unit normal vector of the boundary surface; and ${\mathrm{\theta} _{\textrm{i}}}$ and ${{{\textrm{N}_{\textrm{i}}} \equiv {\mathrm{{\mathrm{\xi}}} _{\textrm{i-1}}}} \mathord{\left/ {\vphantom {{{\textrm{N}_{\textrm{i}}} \equiv {\mathrm{{\mathrm{\xi}}}_{\textrm{i-1}}}} {{\mathrm{{\mathrm{\xi}}}_{\textrm{i}}}}}} \right.} {{\mathrm{{\mathrm{\xi}}} _{\textrm{i}}}}}$ are, respectively, the incidence angle and relative refractive index of medium i with respect to medium i-1.It is noted from Eqs. (5), (8a), and (8b) that the reflected/refracted ray ${\bar{\textrm{R}}_{\textrm{i}}} = {\left[ {\begin{array}{cc} {{{\bar{\textrm{P}}}_{\textrm{i}}}}&{{{\bar{\ell }}_{\textrm{i}}}} \end{array}} \right]^{\textrm{T}}}$ for the given ith boundary surface ${\bar{\textrm{r}}_{\textrm{i}}}$ is a function of the incoming ray ${\overline{\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}}}$. In other words, ${\overline{\textrm{R}}}_{\textrm{i}}$ is a recursive function which calls itself during execution, to give
(9)$${\bar{\textrm{R}}_{\textrm{i}}} = {\bar{\textrm{R}}_{\textrm{i}}}({{{\bar{\textrm{R}}}_{\textrm{i-1}}}} ).$$
3. Taylor series expansion of ${\bar{\textrm{R}}_{\textrm{n}}}$
By successively applying Eq. (9) from i = n to i=1, and then using Eq. (4), the following composite function for ray ${\bar{\textrm{R}}_{\textrm{n}}} = {\left[ {\begin{array}{cc} {{{\bar{\textrm{P}}}_{\textrm{n}}}}&{{{\bar{\ell }}_{\textrm{n}}}} \end{array}} \right]^{\textrm{T}}}$ on the image plane for an optical system with n boundary surfaces is obtained:
(10)$${\bar{\textrm{R}}_{\textrm{n}}}({\bar{\textrm{X}}_0}) = {\bar{\textrm{R}}_{\textrm{n}}}({{{\bar{\textrm{R}}}_{\textrm{n} - \textrm{1}}}({{{\bar{\textrm{R}}}_{\textrm{n} - 2}}({\ldots ({{{\bar{\textrm{R}}}_{\textrm{i}}}({{{\bar{\textrm{R}}}_{\textrm{i-1}}}({\ldots ({{{\bar{\textrm{R}}}_2}({{{\bar{\textrm{R}}}_1}({{{\bar{\textrm{R}}}_0}({{\bar{\textrm{X}}}_0})} )} )} )\ldots } )} )} )} )} )} ).$$
It is noted that the ray ${\bar{\textrm{R}}_{\textrm{n}}}$ is a composite function. A composite function is a function which operates in turn on another function. To evaluate such a function, it is first necessary to evaluate the internal functions, and then determine the outer function based on the results of these internal functions. However, the complete explicit function of Eq. (10) even for a system with a singlet and image plane (which is the simplest system having n=3) is still lengthy, making it is difficult for use in design. Furthermore, it is every difficult to investigate the performance of a system possessing a large number of boundary surfaces using a ray tracing method. For example, for a system with n=40 boundary surfaces, the ray ${\bar{\textrm{R}}_{\textrm{n}}}$ on the image plane has the form of a 40-layered deep composite function. The most elegant solution in mathematics for this type of hard problem is the Taylor series expansion since it can provide the polynomial expansion of ${\bar{\textrm{R}}_{\textrm{n}}}$ in terms of ${\bar{\textrm{X}}_0}$, which is far easier to work with than the function in Eq. (10), i.e., (11)$$\begin{aligned} {{\bar{\textrm{R}}}_{\textrm{n}}}({{\bar{\textrm{X}}}_0}) &= {{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0}) + \frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{X}}}_0}}}{{\bar{\textrm{X}}}_0} + \frac{1}{2}\left[ {{{\left( {{{\bar{\textrm{X}}}_0}} \right)}^{\textrm{T}}}\frac{{{{\mathrm{\partial}} ^2}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} \bar{\textrm{X}}_0^2}} {{\bar{\rm X}}_0}} \right]\\ + &\frac{1}{6}\left[ {{{\left( {{{\bar{\textrm{X}}}_0}} \right)}^{\textrm{T}}}\frac{{{{\mathrm{\partial}} ^3}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} \bar{\textrm{X}}_0^3}} {{\bar{\textrm{X}}}_0}} \right]{{\bar{\textrm{X}}}_0} + \frac{1}{{24}}{\left( {{{\bar{\textrm{X}}}_0}} \right)^{\textrm{T}}}\left[ {{{\left( {{{\bar{\textrm{X}}}_0}} \right)}^{\textrm{T}}}\frac{{{{\mathrm{\partial}} ^4}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} \bar{\textrm{X}}_0^4}} {{\bar{\textrm{X}}}_0}} \right]{{\bar{\textrm{X}}}_0}\\ + &\frac{1}{{120}}\left\{ {{{\left( {{{\bar{\textrm{X}}}_0}} \right)}^{\textrm{T}}}\left[ {{{\left( {{{\bar{\textrm{X}}}_0}} \right)}^{\textrm{T}}}\frac{{{{\mathrm{\partial}} ^5}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} \bar{\textrm{X}}_0^5}} {{\bar{\textrm{X}}}_0}} \right]{{\bar{\textrm{X}}}_0}} \right\}{{\bar{\textrm{X}}}_0} + ....\\ \equiv {{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0}) &+ \Delta {{\bar{\textrm{R}}}_{\textrm{n/1st}}} + \Delta {{\bar{\textrm{R}}}_{\textrm{n/2nd}}} + \Delta {{\bar{\textrm{R}}}_{\textrm{n/3rd}}} + \Delta {{\bar{\textrm{R}}}_{\textrm{n/4th}}} + \Delta {{\bar{\textrm{R}}}_{\textrm{n/5th}}} + ..., \end{aligned}$$
where ${{{{\mathrm{\partial}}^{\textrm{m}}}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})} \mathord{\left/ {\vphantom {{{{\mathrm{\partial}} ^{\textrm{m}}}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})} {{\mathrm{\partial}} \bar{\textrm{X}}_0^{\textrm{m}}}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^{\textrm{m}}}}$ denotes the mth-order derivative of ${\bar{\textrm{R}}_{\textrm{n}}}$ with respect to ${\bar{\textrm{X}}_0}$ and is evaluated at ${\bar{\textrm{X}}_0} = \bar{0}$ (i.e., the reference ray travelling along ${\textrm{z}_0}$ axis, denoted as reference source ray hereafter). The first term, ${\bar{\textrm{R}}_{\textrm{n}}}(\bar{0})$, of Eq. (11) is the intersection point of the image plane and the ray originated from the reference source ray. Meanwhile, the second term, $\Delta {\bar{\textrm{R}}_{\textrm{n/1st}}}$, provides the paraxial ray tracing equations required to predict the image point ${\bar{\textrm{P}}_{\textrm{n}}}$ for a given object point ${\bar{\textrm{P}}_0}$. The remaining terms in Eq. (11) are all associated with the ray aberration functions. Equation (11) is still valid for non-axially symmetrical systems (e.g., system with freeform optical surfaces [14,15]).4. First-order expansion and paraxial ray tracing equations
The first-order expansion of the incidence ray ${\bar{\textrm{R}}_{\textrm{n}}}$ on the image plane with respect to ${\bar{\textrm{X}}_0}$ [i.e., the second term of Eq. (11)] can be determined using the chain rule as
(12)$$\begin{aligned} &\Delta {\bar{\textrm{R}}}_{\textrm{n/1st}} = \frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{X}}}_\textrm{0}}}}{{\bar{\textrm{X}}}_0}\\ & = \frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n - 1}}}}}\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n - 1}}}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{n - 2}}}}}...\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i-1}}}}}...\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_2}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_1}}}\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_1}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_0}}}\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_0}(\bar{0})}}{{{\mathrm{\partial}} {{\bar{\textrm{X}}}_\textrm{0}}}}{{\bar{\textrm{X}}}_0}. \end{aligned}$$
It is seen that Eq. (12) comes from the following equation by setting i = n: (13)$$\begin{aligned} &\Delta {{\bar{\textrm{R}}}_{\textrm{i}}} = \frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i-1}}}}}\Delta {{\bar{\textrm{R}}}_{\textrm{i-1}}}\\ &= \frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i-1}}}}}\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i-1}}}}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i - 2}}}}}\ldots \frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_2}}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_1}}}\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_1}}}{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_0}}}\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_0}}}{{{\mathrm{\partial}} {{\bar{\textrm{X}}}_{0}}}}{{\bar{\textrm{X}}}_0}. \end{aligned}$$
The term ${{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}} {{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i-1}}}}}} \right.} {{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i-1}}}}}$ in Eq. (13) represents the first-order derivative of ray ${\bar{\textrm{R}}_{\textrm{i}}}$ with respect to the incoming ray ${\bar{\textrm{R}}_{\textrm{i-1}}}$. Since the higher-order derivatives of ${\bar{\textrm{R}}_{\textrm{i}}}$ with respect to ${\bar{\textrm{X}}_{0}}$ are bulky, it is inconvenient to use such a notation to express ${{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_{0}}}}} \right.} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_{0}}}}$. Consequently, the following summation notation is used instead: (14)$$\frac{{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}}}{{{\mathrm{\partial}} {{\bar{\textrm{X}}}_{0}}}} = {\left[ {\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right]_{6 \times 5}} = \left[ {\sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right].$$
To better explain the notation of Eq. (14), Fig. 2 shows the structure of the first-order derivative of a vector function $\bar{\textrm{F}}(\bar{\textrm{X}}) = {\left[ {\begin{array}{cc} {{\textrm{f}_{1}}(\bar{\textrm{X}})}&{{\textrm{f}_{2}}(\bar{\textrm{X}})} \end{array}} \right]^{\textrm{T}}} \equiv [{\textrm{f}_{\textrm{q}}} ]$ with $\bar{\textrm{X}} = {\left[ {\begin{array}{ccc} {{\textrm{x}_1}}&{{\textrm{x}_2}}&{{\textrm{x}_3}} \end{array}} \right]^{\textrm{T}}} \equiv [{\textrm{x}_{\textrm{u}}} ]$. Note that the post-subscript “u” in Eq. (14) indicates that ${{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right.} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}$ (i.e., the (q,u)th component of ${{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {{\bar{\textrm{R}}}_{\textrm{i}}}} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_{0}}}}} \right.} {{\mathrm{\partial}} {{\bar{\textrm{X}}}_{0}}}}$) is the derivative of the qth component of ${\bar{\textrm{R}}_{\textrm{i}}}$ with respect to the uth component of ${\bar{\textrm{X}}_{0}}$.The applications of Eq. (13) include determining the image orientation in prism systems (chap. 9 of [3]); and the caustic surfaces, point spread function, and modulation transfer function (chap. 13 of [3]) in non-axially symmetrical optical systems. However, one of its most important applications lies in the paraxial optics field.
In paraxial optics, symbols without and with a prime, respectively, are used to represent the parameters before and after a boundary surface (or system), respectively. The most common examples are ${\textrm{h}_0}$ and ${\textrm{h}^{\prime}_{0}}$, which denote the object height and image height, respectively; or ${\bar{\textrm{Q}_{\textrm{i}}}}$ and ${\bar{\textrm{Q}}^{\prime}_{\textrm{i}}}$, which designate the two points of a ray path before and after the ith paraxial surface ${\bar{\textrm{r}}_{\textrm{i}}}$. In Fig. 3 a paraxial skew ray ${\left[ {{{\begin{array}{cc} {{{\bar{\textrm{P}}^{\prime}}_{\textrm{i-1}}}}&{\bar{\ell }^{\prime}} \end{array}}_{\textrm{i-1}}}} \right]^{\textrm{T}}}$ originates at ${\bar{\textrm{Q}^{\prime}_{\textrm{i-1}}}}$ and propagates along a straight-line path in a medium with refractive index ${\upsilon _{\textrm{i}}}$ until it reaches point ${\bar{\textrm{Q}}_{\textrm{i}}}$. The term $\Delta {\bar{\textrm{R}}_{\textrm{i}}} = {\left[ {\begin{array}{cc} {\Delta {{\bar{\textrm{P}}}_{\textrm{i}}}}&{\Delta {{\bar{\ell }}_{\textrm{i}}}} \end{array}} \right]^{\textrm{T}}}$ of Eq. (13) is the differential ray between a skew paraxial ray ${\bar{\textrm{R}}_{\textrm{i}}} = {\left[ {\begin{array}{cc} {{{\bar{\textrm{P}}}_{\textrm{i}}}}&{{{\bar{\ell }}_{\textrm{i}}}} \end{array}} \right]^{\textrm{T}}}$ at ${\bar{\textrm{Q}}_{\textrm{i}}}$ and the reference ray ${\bar{\textrm{R}}_{\textrm{i,axis}}} = {\left[ {\begin{array}{cc} {{{\bar{\textrm{P}}}_{\textrm{i,axis}}}}&{{{\bar{\ell }}_{\textrm{i,axis}}}} \end{array}} \right]^{\textrm{T}}}$ originated from reference source ray. It then undergoes refraction on ${\bar{\textrm{r}}_{\textrm{i}}}$ and propagates again through the new medium with refractive index ${\upsilon ^{\prime}_{\textrm{i}}}$ to a neighboring point ${\bar{\textrm{Q}}^{\prime}_{\textrm{i}}}$ after ${\bar{\textrm{r}}_{\textrm{i}}}$. The differential ray, ${\left[ {{{\begin{array}{cc} {\Delta {{\bar{\textrm{P}}^{\prime}}_{\textrm{i}}}}&{\Delta \bar{\ell }^{\prime}} \end{array}}_{\textrm{i}}}} \right]^{\textrm{T}}}$, at ${\bar{\textrm{Q}^{\prime}_{\textrm{i}}}}$ can be estimated via the following $6 \times 6$ matrix manipulation (Eq. (4.12) of [3]):
(15)$$\left[ {\begin{array}{c} {\Delta {{\bar{\textrm{P}}^{\prime}_{\textrm{i}}}}}\\ {\Delta {{\bar{\ell }^{\prime}_{\textrm{i}}}}} \end{array}} \right] = {\bar{\textrm{M}}_{\textrm{i}}}{\bar{\textrm{T}}_{\textrm{i}}}\left[ {\begin{array}{c} {\Delta {{\bar{\textrm{P}}^{\prime}_{\textrm{i-1}}}}}\\ {\Delta {{\bar{\ell }^{\prime}_{\textrm{i-1}}}}} \end{array}} \right],$$
where ${\bar{\textrm{M}}_{\textrm{i}}}$ and ${\bar{\textrm{T}}_{\textrm{i}}}$ are the transfer and refraction matrices and are given respectively as (16a)$${\bar{\textrm{M}}_{\textrm{i}}} = \left[ \begin{array}{cccccc} 1 &0 &0 &0 &0 &0\\ 0 &0 &0 &0 &0 &0\\ 0 &0 &1 &0 &0 &0\\ {{{({\textrm{N}_{\textrm{i}}} - 1)} \mathord{\left/ {\vphantom {{({\textrm{N}_{\textrm{i}}} - 1)} {{\textrm{R}_{\textrm{i}}}}}} \right.} {{\textrm{R}_{\textrm{i}}}}}} &0 &0 &{{\textrm{N}_{\textrm{i}}}} &0 &0\\ 0 &0 &0 &0 &{\textrm{N}_{\textrm{i}}^2} &0\\ 0 &0 &{{{({\textrm{N}_{\textrm{i}}} - 1)} \mathord{\left/ {\vphantom {{({\textrm{N}_{\textrm{i}}} - 1)} {{\textrm{R}_{\textrm{i}}}}}} \right.} {{\textrm{R}_{\textrm{i}}}}}} &0 &0 &{{\textrm{N}_{\textrm{i}}}} \end{array} \right],$$
(16b)$${\bar{\textrm{T}}_{\textrm{i}}} = \left[ \begin{array}{cccccc} 1 &0 &0 &{{\lambda_{\textrm{i,axis}}}} &0 &0\\ 0 &1 &0 &0 &{{\lambda_{\textrm{i,axis}}}} &0\\ 0 &0 &1 &0 &0 &{{\lambda_{\textrm{i,axis}}}}\\ 0 &0 &0 &1 &0 &0\\ 0 &0 &0 &0 &1 &0\\ 0 &0 &0 &0 &0 &1 \end{array} \right].$$
Equation (15) should be performed in a local coordinate frame ${(\textrm{xyz})_{\textrm{i}}}$ aligned such that the ${\textrm{z}_{\textrm{i}}}$ axis coincides with the optical axis with length ${\lambda _{\textrm{i,axis}}}$. Term ${\textrm{N}_{\textrm{i}}} = {{{\upsilon _{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\upsilon_{\textrm{i}}}} {{{\upsilon^{\prime}_{\textrm{i}}}}}}} \right.} {{{\upsilon ^{\prime}_{\textrm{i}}}}}}$ in Eq. (16a) is the relative refraction index. Equation (16a) is valid for both convex (i.e., positive value of ${\textrm{R}_{\textrm{i}}}$) and concave (i.e., negative value of ${\textrm{R}_{\textrm{i}}}$) boundary surfaces, as well as flat boundaries (${\textrm{R}_{\textrm{i}}} = \infty $).
For an axis-symmetrical system, the object point ${\bar{\textrm{P}}_0}$ can be positioned without any loss of generality on a meridional plane located at ${\textrm{P}_{0{\textrm{z}}}}$ with a height ${\textrm{h}_0}$ by setting ${\textrm{P}_{\textrm{0x}}} = {\mathrm{\alpha} _0} = 0$ in Eq. (3). To perform the paraxial ray tracing for an axis-symmetrical system, it is first necessary to use the height-slope matrices, ${\left[ {\begin{array}{cc} {{\textrm{y}_{\textrm{i}}}}&{{\upsilon_{\textrm{i}}}{\mathrm{\mu}_{\textrm{i}}}} \end{array}} \right]^{\textrm{T}}}$ and ${\left[ {\begin{array}{cc} {{{{\rm y}^{\prime}_{\textrm{i}}}}}&{{{\upsilon^{\prime}_{\textrm{i}}}}{{\mathrm{\mu}^{\prime}_{\textrm{i}}}}} \end{array}} \right]^{\textrm{T}}}$, to identify points ${\bar{\textrm{Q}}_{\textrm{i}}}$ and ${\bar{\textrm{Q}}^{\prime}_{\textrm{i}}}$, respectively (Fig. 4). Note that ${\textrm{y}_{\textrm{i}}}$ and ${\mathrm{\mu} _{\textrm{i}}}$ respectively denote the height of ${\bar{\textrm{Q}}_{\textrm{i}}}$ and the paraxial angle ${\mathrm{\mu} _{\textrm{i}}}$ the ray makes with the optical axis. For an axis-symmetrical system, Eqs. (15), (16a), and (16b) can be simplified as Eqs. (4.18) and (4.19) of [3].
(17a)$$\left[ {\begin{array}{c} {{\textrm{y}_{\textrm{i}}}}\\ {{\upsilon_{\textrm{i}}}{\mathrm{\mu}_{\textrm{i}}}} \end{array}} \right] = {\bar{\textrm{M}}_{\textrm{i}}}{\bar{\textrm{T}}_{\textrm{i}}}\left[ {\begin{array}{c} {{\textrm{y}_{\textrm{i-1}}}}\\ {{\upsilon_{\textrm{i-1}}}{\mathrm{\mu}_{\textrm{i-1}}}} \end{array}} \right],$$
(17b)$${\bar{\textrm{T}}_{\textrm{i}}} = \left[ {\begin{array}{cc} 1 &{{{{\lambda_{\textrm{i,axis}}}} \mathord{\left/ {\vphantom {{{\lambda_{\textrm{i,axis}}}} {{\upsilon_{\textrm{i}}}}}} \right.} {{\upsilon_{\textrm{i}}}}}}\\ 0 &1 \end{array}} \right],$$
(17c)$${\bar{\textrm{M}}_{\textrm{i}}} = \left[ {\begin{array}{cc} 1 &0\\ {{{({\upsilon_{\textrm{i}}} - {{\upsilon^{\prime}_{\textrm{i}}}})} \mathord{\left/ {\vphantom {{({\upsilon_{\textrm{i}}} - {{\upsilon^{\prime}_{\textrm{i}}}})} {{\textrm{R}_{\textrm{i}}}}}} \right.} {{\textrm{R}_{\textrm{i}}}}}} &1 \end{array}} \right].$$
Equation (17a) describes the path of a paraxial meridional ray traveling along the ${\textrm{y}_0}{\textrm{z}_0}$ plane through an axis-symmetrical system. It can also be used to determine the six cardinal points, the image position, and the lateral and longitudinal magnifications of an object close to the optical axis. Furthermore, it is possible by using Eq. (17a) to trace the marginal and chief paraxial rays to have Buchdahl aberration coefficients [9].
5. Higher-order expansions and ray aberrations
Further differentiation of Eq. (14) with respect to ${\bar{\textrm{X}}_{0}}$ yields the following second-order derivative matrix of ray ${\bar{\textrm{R}}_{\textrm{i}}}$ with respect to ${\bar{\textrm{X}}_{0}}$:
(18)$$\begin{aligned} \frac{{{\mathrm{\partial}} {}^2{{\bar{\textrm{R}}}_{\textrm{i}}}}}{{{\mathrm{\partial}} \bar{\textrm{X}}_0^2}} = &{\left[ {\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right]_{6 \times 5 \times 5}}\\ = &\left[ {\sum\limits_{\textrm{g} = 1}^6 {\left( {\sum\limits_{\textrm{f} = 1}^6 {\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} } \right)} + \sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right], \end{aligned}$$
where the post-subscript “v” is the second-order derivative index. To better explain the notation of Eq. (18), Fig. 5 illustrates the architecture of the second-order derivative of a vector function $\bar{\textrm{F}}(\bar{\textrm{X}},\bar{\textrm{Y}}\textrm{) = }{\left[ {\begin{array}{cc} {{\textrm{f}_{1}}(\bar{\textrm{X}},\bar{\textrm{Y}})}&{{\textrm{f}_{2}}(\bar{\textrm{X}},\bar{\textrm{Y}})} \end{array}} \right]^{\textrm{T}}} \equiv [{\textrm{f}_{\textrm{u}}} ]$ with $\bar{\textrm{X}} = {\left[ {\begin{array}{ccc} {{\textrm{x}_1}}&{{\textrm{x}_2}}&{{\textrm{x}_3}} \end{array}} \right]^{\textrm{T}}} \equiv [{\textrm{x}_{\textrm{u}}}]$ and $\bar{\textrm{Y}} = {\left[ {\begin{array}{cc} {{\textrm{y}_1}}&{{\textrm{y}_2}} \end{array}} \right]^{\textrm{T}}} \equiv [{\textrm{y}_{\textrm{v}}} ]$. Note that the required first-order derivatives of Eq. (18) (i.e., ${{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}}}} \right.} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}}}$, ${{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right.} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}$ and ${{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}} {{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}}}} \right.} {{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}$) are obtained using Eq. (14). The final term of Eq. (18) (i.e., ${{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{i-1f}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{i-1f}}}} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right.} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}$) is known since we have ${{{\mathrm{\partial}} {}^2{{\bar{\textrm{R}}}_0}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^2{{\bar{\textrm{R}}}_0}} {{\mathrm{\partial}} \bar{\textrm{X}}_0^2}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^2}}$ from Eqs. (1) and (2). Equation (18) is then computed sequentially from i=1 to i = n to give ${{{{\mathrm{\partial}} ^2}{{\bar{\textrm{R}}}_{\textrm{n}}}} \mathord{\left/ {\vphantom {{{{\mathrm{\partial}}^2}{{\bar{\textrm{R}}}_{\textrm{n}}}} {{\mathrm{\partial}} \bar{\textrm{X}}_0^2}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^2}}$. The numerical results show that ${{{{\mathrm{\partial}} ^2}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})} \mathord{\left/ {\vphantom {{{{\mathrm{\partial}}^2}{{\bar{\textrm{R}}}_{\textrm{n}}}(\bar{0})} {{\mathrm{\partial}} \bar{\textrm{X}}_0^2}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^2}}$ is a zero vector for an axis-symmetrical system. Consequently, $\Delta {\bar{\textrm{R}}_{\textrm{n/2nd}}}$ in Eq. (11) provides no contribution to its ray aberrations. However, it may introduce aberration in a non-axially symmetrical system. One example of such systems is the imaging system with freeform surfaces.Further differentiating Eq. (18) with respect to ${\bar{\textrm{X}}_{0}}$ gives the third-order derivative matrix of ray ${\bar{\textrm{R}}_{\textrm{i}}}$ with respect to ${\bar{\textrm{X}}_{0}}$, i.e.,
(19)$$\begin{aligned} \frac{{{\mathrm{\partial}} {}^3{{\bar{\textrm{R}}}_{\textrm{i}}}}}{{{\mathrm{\partial}} \bar{\textrm{X}}_0^3}} &= {\left[ {\frac{{{\mathrm{\partial}} {}^3{\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right]_{6 \times 5 \times 5 \times 5}}\\ = &\left[ {\sum\limits_{\textrm{h} = 1}^6 {\left( {\sum\limits_{\textrm{g} = 1}^6 {\left( {\sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {}^3{\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1h}}}{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1h}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right)} } \right)} } \right.\\ &+ \sum\limits_{\textrm{g} = 1}^6 {\left( {\sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{i-1g}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right)} \\ &+ \sum\limits_{\textrm{g} = 1}^6 {\left( {\sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}}}\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right)} \\ &\left. { + \sum\limits_{\textrm{g} = 1}^6 {\left( {\sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1g}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}}}\frac{{{\mathrm{\partial}} {}^2{\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right)} + \sum\limits_{\textrm{f} = 1}^6 {\left( {\frac{{{\mathrm{\partial}} {\textrm{R}_{\textrm{iq}}}}}{{{\mathrm{\partial}} {\textrm{R}_{\textrm{i-1f}}}}}\frac{{{\mathrm{\partial}} {}^3{\textrm{R}_{\textrm{i-1f}}}}}{{{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right)} } \right], \end{aligned}$$
where the post-subscript “w” represents the third-order derivative index. Note that the required first- and second-order derivatives of Eq. (19) are extracted from Eqs. (14) and (18), respectively. The last term (i.e., ${{{\mathrm{\partial}} {}^3{\textrm{R}_{\textrm{i-1f}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^3{\textrm{R}_{\textrm{i-1f}}}} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}} \right.} {{\mathrm{\partial}} {\textrm{X}_{\textrm{0w}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0v}}}{\mathrm{\partial}} {\textrm{X}_{\textrm{0u}}}}}$) of Eq. (19) is known since ${{{\mathrm{\partial}} {}^3{{\bar{\textrm{R}}}_0}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^3{{\bar{\textrm{R}}}_0}} {{\mathrm{\partial}} \bar{\textrm{X}}_0^3}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^3}}$ is obtained from Eqs. (1) and (2). Equation (19) is then determined sequentially from i=1 to i = n to obtain the term $\Delta {\bar{\textrm{R}}_{\textrm{n/3rd}}}$ in Eq. (11). The most attractive application of Eq. (19) is the determination of Seidel primary ray aberration coefficients.To have the Seidel primary ray aberration coefficients for an axis-symmetrical system, the object point ${\overline{\textrm{P}}_0} = {\left[ {\begin{array}{ccc} 0 &{\textrm{h}_0} &{\textrm{P}_{0{\textrm{z}}}} \end{array}} \right]^{\textrm{T}}}$ can be placed on the meridional plane without any loss of generality in order to obtain the Seidel variable vector ${\left[ {\begin{array}{ccccc} 0 &{{\textrm{h}_0}} &{{\textrm{P}_{\textrm{0z}}}} &{{\mathrm{\alpha}_0}} &{{\mathrm{\beta}_0}} \end{array}} \right]^{\textrm{T}}}$. It is noted that only the first and second components of $\Delta {\bar{\textrm{R}}_{\textrm{n/3rd}}}$, $\Delta {\textrm{P}_{\textrm{nx/3rd}}}$ and $\Delta {\textrm{P}_{\textrm{ny/3rd}}}$, are required to investigate the primary Seidel ray aberration coefficients. Due to the axis-symmetrical nature of the considered system, various terms in ${{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}(\bar{0})} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}(\bar{0})} {{\mathrm{\partial}} \bar{\textrm{X}}_0^3}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^3}}$ and ${{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}(\bar{0})} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}(\bar{0})} {{\mathrm{\partial}} \bar{\textrm{X}}_0^3}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^3}}$ disappear, to leave
(20a)$$\begin{aligned} \Delta {\textrm{P}_{\textrm{nx/3rd}}} &= \frac{1}{6}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha} _0^3}}\mathrm{\alpha} _0^3 + \frac{1}{2}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} {\mathrm{\alpha} _0}{\mathrm{\partial}} \mathrm{\beta}_0^2}}{\mathrm{\alpha} _0}\mathrm{\beta} _0^2\\ + &\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} {\textrm{h}_0}{\mathrm{\partial}} {\mathrm{\alpha}_0}{\mathrm{\partial}} {\mathrm{\beta} _0}}}{\textrm{h}_0}{\mathrm{\alpha} _0}{\mathrm{\beta} _0} + \frac{1}{2}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} {\textrm{h}}_0^2{\mathrm{\partial}} {\mathrm{\alpha} _0}}}\textrm{h}_0^2{\mathrm{\alpha} _0}, \end{aligned}$$
(20b)$$\begin{aligned} \Delta {\textrm{P}_{\textrm{ny/3rd}}} &= \frac{1}{6}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}}}{{{\mathrm{\partial}} \mathrm{\beta} _0^3}}\mathrm{\beta} _0^3 + \frac{1}{2}\frac{{{\mathrm{\partial}} {}^3{P_{\textrm{ny}}}}}{{{\mathrm{\partial}} \mathrm{\alpha} _0^2{\mathrm{\partial}} {\mathrm{\beta} _0}}}\mathrm{\alpha}_0^2{\mathrm{\beta} _0} + \frac{1}{2}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}}}{{{\mathrm{\partial}} {\textrm{h}_0}{\mathrm{\partial}} \mathrm{\alpha} _0^2}}{\textrm{h}_0}\mathrm{\alpha} _0^2\\ + &\frac{1}{2}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}}}{{{\mathrm{\partial}} {\textrm{h}_0}{\mathrm{\partial}} \mathrm{\beta} _0^2}}{\textrm{h}_0}\mathrm{\beta} _0^2 + \frac{1}{2}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}}}{{{\mathrm{\partial}} \textrm{h}_0^2{\mathrm{\partial}} {\mathrm{\beta} _0}}}\textrm{h}_0^2{\mathrm{\beta} _0} + \frac{1}{6}\frac{{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{ny}}}}}{{{\mathrm{\partial}} \textrm{h}_0^3}}\textrm{h}_0^3.\end{aligned}$$
The coefficients of spherical aberration, coma, astigmatism, field curvature and distortion (denoted as ${\textrm{B}_{1}}$, ${\textrm{B}_{2}}$, ${\textrm{B}_{3}}$, ${\textrm{B}_{4}}$ and ${\textrm{B}_{5}}$, respectively) can be obtained from Eqs. (20a) and (20b) (see Eqs. (22)–(28) of [13]). The numerical results of [13] show that the five aberration coefficients thus determined are in excellent agreement with those obtained by Zemax simulations.It is possible to have higher order derivative matrix of ${\bar{\textrm{R}}_{\textrm{i}}}$ [e.g., the fourth- and fifth-order derivative matrices, ${{{\mathrm{\partial}} {}^4{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^4{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\mathrm{\partial}} \bar{\textrm{X}}_0^4}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^4}}$ and ${{{\mathrm{\partial}} {}^5{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^5{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\mathrm{\partial}} \bar{\textrm{X}}_0^5}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^5}}$, in Eq. (11)] by further differentiating Eq. (19) with respect to ${\bar{\textrm{X}}_{0}}$. However, determining these higher order ray derivative matrices is extremely challenging. Note that, due to the axis-symmetrical nature of an axis-symmetrical system, the fourth-order derivative matrix ${{{\mathrm{\partial}} {}^4{{\bar{\textrm{R}}}_{\textrm{i}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^4{{\bar{\textrm{R}}}_{\textrm{i}}}} {{\mathrm{\partial}} \bar{\textrm{X}}_0^4}}} \right.} {{\mathrm{\partial}} \bar{\textrm{X}}_0^4}}$ provides no contribution to ray aberrations. To investigate the secondary ray aberration for such system, its fifth-order of the ray derivative (e.g., ${{{{\mathrm{\partial}} ^5}{\textrm{P}_{\textrm{nx}}}} \mathord{\left/ {\vphantom {{{{\mathrm{\partial}}^5}{\textrm{P}_{\textrm{nx}}}} {{\mathrm{\partial}} \mathrm{\alpha}_0^5}}} \right.} {{\mathrm{\partial}} \mathrm{\alpha} _0^5}}$) can be approximated as
(21)$$\frac{{{{\mathrm{\partial}} ^5}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha} _0^5}} \approx \frac{{ - {{\left( {\frac{{{{\mathrm{\partial}}^3}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha}_0^3}}} \right)}_2} + 16{{\left( {\frac{{{{\mathrm{\partial}}^3}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha}_0^3}}} \right)}_1} - 30{{\left( {\frac{{{{\mathrm{\partial}}^3}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha}_0^3}}} \right)}_0} + 16{{\left( {\frac{{{{\mathrm{\partial}}^3}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha}_0^3}}} \right)}_{ - 1}} - {{\left( {\frac{{{{\mathrm{\partial}}^3}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha}_0^3}}} \right)}_{ - 2}}}}{{12{{(\Delta {\mathrm{\alpha} _0})}^2}}},$$
where the fifth-order derivatives ${{{{\mathrm{\partial}} ^5}{\textrm{P}_{\textrm{nx}}}} \mathord{\left/ {\vphantom {{{{\mathrm{\partial}}^5}{\textrm{P}_{\textrm{nx}}}} {{\mathrm{\partial}} \mathrm{\alpha}_0^5}}} \right.} {{\mathrm{\partial}} \mathrm{\alpha} _0^5}}$ is estimated from the results of the five third-order derivatives (i.e., ${{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}} \mathord{\left/ {\vphantom {{{\mathrm{\partial}} {}^3{\textrm{P}_{\textrm{nx}}}} {{\mathrm{\partial}} \mathrm{\alpha}_0^3}}} \right.} {{\mathrm{\partial}} \mathrm{\alpha} _0^3}}$) evaluated by five neighboring rays of optical axis. Our study shows that the ${\textrm{C}_{1}}$ coefficient (Eq. (3.2) in p. 62 of [1]) of the secondary spherical aberration can be determined by (22)$${\textrm{C}_{1}} = \frac{1}{{120}}\frac{{{{\mathrm{\partial}} ^5}{\textrm{P}_{\textrm{nx}}}}}{{{\mathrm{\partial}} \mathrm{\alpha} _0^5}}{(\frac{{{\mathrm{\alpha} _0}}}{{{\textrm{x}_{\textrm{a}}}}})^5} = - 2.1865274 \times \textrm{1}{\textrm{0}^{ - 8}},$$
where $\textrm{x}_{\textrm{a}}$ is the coordinate of the entrance pupil (see Fig. 2 of [13]), and ${{{{\mathrm{\partial}} ^5}{\textrm{P}_{\textrm{nx}}}} \mathord{\left/ {\vphantom {{{{\mathrm{\partial}}^5}{\textrm{P}_{\textrm{nx}}}} {{\mathrm{\partial}} \mathrm{\alpha}_0^5}}} \right.} {{\mathrm{\partial}} \mathrm{\alpha} _0^5}} ={-} 1418113.3208$, ${{{\mathrm{\alpha} _0}} \mathord{\left/ {\vphantom {{{\mathrm{\alpha}_0}} {{\textrm{x}_{\textrm{a}}} = 4.502410 \times {{10}^{ - 3}}}}} \right.} {{\textrm{x}_{\textrm{a}}} = 4.502410 \times {{10}^{ - 3}}}}$. The spherical aberration coefficient is ${\textrm{C}_1}\mathrm{\rho} _{\max }^{5} = - 0.089300$ for the system of Fig. 3 in [13] (with an entrance pupil having maximum opening radius ${\mathrm{\rho} _{\max }} = 21$ and the on-axis object is at ${\textrm{P}_{\textrm{0z}}} ={-} 200$). This value has $15.3\%$ percentage error compared with the value $- 0.075594$ from Zemax.6. Conclusions
This paper has proposed a roadmap for geometrical optics based on the Taylor series expansion of a skew ray with respect to the source ray variable vector ${\bar{\textrm{X}}_0}$. It has been indicated that the paraxial ray tracing equations, point spread function, caustic surfaces, modulation transfer functions, and prism image orientation can all be explored from its first-order expansion, while the primary ray aberrations can be determined from the third-order expansion. We have computed the Seidel primary ray aberrations, and appears to offer the opportunity to compute the higher-order ray aberrations simply by generating higher-order ray derivative matrices. Since the derivatives of the optical path length (OPL) are supplementary quantities of the ray derivatives, the proposed method also provides a feasible approach for determining the wavefront aberrations by using the Taylor series expansion of OPL, which can be obtained by replacing ${\bar{\textrm{R}}_{\textrm{n}}}$ of Eq. (11) with OPL. Furthermore, it is possible to determine the wavefront aberrations based on Zernike polynomials, provided polar coordinate systems are used. Thus, it is speculated that the Taylor series expansion may serve as a useful basis for the development of new approaches to explore a wide variety of problems in the geometrical optics field.
Funding
Ministry of Science and Technology, Taiwan (106-2221-E-006-091-MY3).
Disclosures
The author declare no conflicts of interest.
References
1. W. J. Smith, Modern Optical Engineering, 3rd ed. (McGraw-Hill, 2000)
2. M. Avendaño-Alejo, D. González-Utrera, and L. Castañeda, “Caustics in a meridional plane produced by plano-convex conic lenses,” J. Opt. Soc. Am. A 28(12), 2619–2628 (2011). [CrossRef]
3. P. D. Lin, Advanced Geometrical Optics (Springer, 2016).
4. V. N. Mahajan, Optical Imaging and Aberrations part I Ray Geometrical Optics (SPIE- The Interational Society for Optical Engineering, 1998).
5. C. S. Liu and P. D. Lin, “computational method for deriving the geometrical point spread function of an optical system,” Appl. Opt. 49(1), 126–136 (2010). [CrossRef]
6. D. G. Burkhard and D. L. Shealy, “Simplified formula for the illuminance in an optical system,” Appl. Opt. 20(5), 897–909 (1981). [CrossRef]
7. P. D. Lin and C. K. Sung, “Matrix-based paraxial skew ray-tracing in 3D systems with non-coplanar optical axis,” Optik 117(7), 329–340 (2006). [CrossRef]
8. R. Kingslake and R. B. Johnson, Lens Design Fundamentals, 2 Edition (Academic, 2010).
9. H. A. Buchdahl, Optical Aberration Coefficients (Dover, 1968).
10. W. T. Welford, Aberrations of Optical Systems (Adam Hilger, 1986).
11. B. Chen and A. M. Herkommer, “High order surface aberrations contributions from phase space analysis of differential rays,” Opt. Express 24(6), 5934–5945 (2016). [CrossRef]
12. J. Sasián, Introduction to Aberrations in Optical Imaging Systems (Cambridge, 2013).
13. P. D. Lin and R. B. Johnson, “Seidel aberration coefficients: an alternative computational method,” Opt. Express 27(14), 19712–19725 (2019). [CrossRef]
14. J. C. Valencia-Estrada and J. García-Márquez, “Freeform geometrical optics I: principles,” Appl. Opt. 58(34), 9455–9464 (2019). [CrossRef]
15. T. Grillon, C. Valencia-Estrada, J. Garcia-Márquez, A. Espinoza-Garcia, and B. Béchadergue, “Freeform geometrical optics II: from parametric representation to CAD/CAM,” Appl. Opt. 58(34), 9465–9472 (2019). [CrossRef]
Open Access
Add to BibSonomy
