Abstract

Differential ray tracing determines an optical system's first-order properties by finding the first-order changes in the configuration of an exiting ray in terms of changes in that ray's initial configuration. When one or more of the elements of a system is inhomogeneous, the only established procedure for carrying out a first-order analysis of a general ray uses relatively inefficient finite differences. To trace a ray through an inhomogeneous medium, one must, in general, numerically integrate an ordinary differential equation, and Runge–Kutta schemes are well suited to this application. We present an extension of standard Runge–Kutta schemes that gives exact derivatives of the numerically approximated rays.

© 1997 Optical Society of America

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  1. A general discussion of differential ray tracing can be found in H. A. Buchdahl, An Introduction to Hamiltonian Optics (Dover, New York, 1993), Sec. 10.
  2. See, for example, R. Kingslake, Lens Design Fundamentals (Academic, Orlando, Fla., 1978), Chap. 10, Sec. I.
  3. See, for example, D. P. Feder, “Differentiation of ray-tracing equations with respect to constructionparameters of rotationally symmetric optics,” J. Opt. Soc. Am. 58, 1494–1505 (1968).
    [CrossRef]
  4. For the mechanics of differential ray tracing through homogeneous media see, for example, A. Cox, A System of Optical Design (Focal Press, London, 1964), pp. 112–121.
  5. E. W. Marchand, Gradient Index Optics (Academic, New York, 1978), Sec. 6.2.
  6. P. J. Sands, “Inhomogeneous lenses. III. Paraxial optics,” J. Opt. Soc. Am. 61, 879–885 (1971).
    [CrossRef]
  7. S. Doric, “Paraxial ray trace for rotationally symmetric homogeneous and inhomogeneousmedia,” J. Opt. Soc. Am. A 1, 818–821 (1984).
    [CrossRef]
  8. K. Tanaka, “Paraxial theory of rotationally distributed-index media by means ofGaussian constants,” Appl. Opt. 23, 1700–1706 (1984).
    [CrossRef]
  9. Other coordinates may be chosen, and the principles presented in Subsection 2.D can still be applied. The use of Cartesian coordinates and optical direction cosines to specify position and direction is convenient for the particular method for numerical integration used here.
  10. More generally, the matrix ∂v/∂w, where v= (v1, v2, ⋯ , vj) and w=(w1, w2, ⋯ , wk), is shorthand for the j×k matrix given by ∂v∂w=∂v1∂w1∂v1∂w2⋯∂v1∂wk∂v2∂w1∂v2∂w2⋮⋮⋱∂vj∂w1⋯∂vj∂wk.
  11. The outer product of two vectors, v=(v1, v2, ⋯ , vj) and w=(w1, w2, ⋯ , wk), is the j×k matrix given by v⊗w=v1w1v1w2⋯v1wkv2w1v2w2⋮⋱⋮vjw1⋯vjwk.
  12. Specifically, J is given by J=[nα′ cos(Iα′)-nα cos(Iα)]1+∂xα∂yα21/2, where Iα is the angle of incidence and Iα′ is the angle of refraction.
  13. See, for example, J. L. Synge, Geometrical Optics—An Introduction to Hamilton's Method, Vol. XXXVII of Cambridge Tracts in Mathematics and Mathematical Physics, F. Smithies, J. A. Todd, eds. (Cambridge U. Press, London, 1962), Chap. 5, Sec. 20.
  14. A. Sharma, D. V. Kumar, A. K. Ghatak, “Tracing rays through graded-index media: a new method,” Appl. Opt. 21, 984–987 (1982).
    [CrossRef] [PubMed]
  15. A Runge–Kutta scheme with Aj defined as in Eq. (3.5) is sometimes referred to as an explicit Runge–Kutta scheme.
  16. B. D. Stone, G. W. Forbes, “Optimal interpolants for Runge–Kutta ray tracing in inhomogeneous media,” J. Opt. Soc. Am. 7, 248–254 (1990).
  17. See the convention given in Note 10 for the meaning of the 3×2 matrices such as ∂Aj/∂y(k).
  18. See, for example, D. M. Young, R. T. Gregory, A Survey of Numerical Mathematics (Dover, New York, 1988), Vol. I, Sec. 7.3.
  19. For a discussion of the use of higher derivatives in numerical computations see R. W. Hamming, Numerical Methods for Scientists and Engineers (Dover, New York, 1986), Chap. 17. Also note that one common way of specifying an index distribution is in the form of a Taylor series. In the most general case, this series can be written as n2(x, y, z)=∑i=0mi∑j=0mj∑k=0mkNijkxiyjzk. In this form, it is less costly to determine higher derivatives of the index than the index itself.
  20. For example, with 18 extra function evaluations required for finite differencing, a simple index function whose evaluation requires five or six operations (if the index is represented in the form of the series given in Note 19, only a couple of terms require this many operations) would make the number of operations for implementing the specific form of finite differencing discussed here and the differential Runge–Kutta scheme approximately equal. With 10 operations required for evaluating the index (or one of its derivatives), the differential Runge–Kutta scheme would require fewer operations.
  21. J. B. Caldwell, D. T. Moore, “Design of gradient-index lens systems for disc format cameras,” Appl. Opt. 25, 3351–3355 (1986), Lens GT2.
    [CrossRef] [PubMed]
  22. B. D. Stone, G. W. Forbes, “First order properties of general optical systems,” J. Opt. Soc. Am. A 9, 478–489 (1992).
    [CrossRef]
  23. Although eight significant figures may seem excessive, they are used so that anyone trying to implement some form of differential ray tracing through inhomogeneous media can use these results to determine whether the program is performing properly.
  24. This follows, for example, from the reference cited in Note 1. Note that determining the error in the determinant of M(α,β) provides a simple check of numerical results.
  25. G. W. Forbes, “Accuracy doubling in the determination of the final ray configurations,” J. Opt. Soc. Am. A 6, 1776–1783 (1989).
    [CrossRef]

1992 (1)

1990 (1)

B. D. Stone, G. W. Forbes, “Optimal interpolants for Runge–Kutta ray tracing in inhomogeneous media,” J. Opt. Soc. Am. 7, 248–254 (1990).

1989 (1)

1986 (1)

1984 (2)

1982 (1)

1971 (1)

1968 (1)

Buchdahl, H. A.

A general discussion of differential ray tracing can be found in H. A. Buchdahl, An Introduction to Hamiltonian Optics (Dover, New York, 1993), Sec. 10.

Caldwell, J. B.

Cox, A.

For the mechanics of differential ray tracing through homogeneous media see, for example, A. Cox, A System of Optical Design (Focal Press, London, 1964), pp. 112–121.

Doric, S.

Feder, D. P.

Forbes, G. W.

Ghatak, A. K.

Gregory, R. T.

See, for example, D. M. Young, R. T. Gregory, A Survey of Numerical Mathematics (Dover, New York, 1988), Vol. I, Sec. 7.3.

Hamming, R. W.

For a discussion of the use of higher derivatives in numerical computations see R. W. Hamming, Numerical Methods for Scientists and Engineers (Dover, New York, 1986), Chap. 17. Also note that one common way of specifying an index distribution is in the form of a Taylor series. In the most general case, this series can be written as n2(x, y, z)=∑i=0mi∑j=0mj∑k=0mkNijkxiyjzk. In this form, it is less costly to determine higher derivatives of the index than the index itself.

Kingslake, R.

See, for example, R. Kingslake, Lens Design Fundamentals (Academic, Orlando, Fla., 1978), Chap. 10, Sec. I.

Kumar, D. V.

Marchand, E. W.

E. W. Marchand, Gradient Index Optics (Academic, New York, 1978), Sec. 6.2.

Moore, D. T.

Sands, P. J.

Sharma, A.

Stone, B. D.

B. D. Stone, G. W. Forbes, “First order properties of general optical systems,” J. Opt. Soc. Am. A 9, 478–489 (1992).
[CrossRef]

B. D. Stone, G. W. Forbes, “Optimal interpolants for Runge–Kutta ray tracing in inhomogeneous media,” J. Opt. Soc. Am. 7, 248–254 (1990).

Synge, J. L.

See, for example, J. L. Synge, Geometrical Optics—An Introduction to Hamilton's Method, Vol. XXXVII of Cambridge Tracts in Mathematics and Mathematical Physics, F. Smithies, J. A. Todd, eds. (Cambridge U. Press, London, 1962), Chap. 5, Sec. 20.

Tanaka, K.

Young, D. M.

See, for example, D. M. Young, R. T. Gregory, A Survey of Numerical Mathematics (Dover, New York, 1988), Vol. I, Sec. 7.3.

Appl. Opt. (3)

J. Opt. Soc. Am. (3)

J. Opt. Soc. Am. A (3)

Other (16)

Although eight significant figures may seem excessive, they are used so that anyone trying to implement some form of differential ray tracing through inhomogeneous media can use these results to determine whether the program is performing properly.

This follows, for example, from the reference cited in Note 1. Note that determining the error in the determinant of M(α,β) provides a simple check of numerical results.

A Runge–Kutta scheme with Aj defined as in Eq. (3.5) is sometimes referred to as an explicit Runge–Kutta scheme.

See the convention given in Note 10 for the meaning of the 3×2 matrices such as ∂Aj/∂y(k).

See, for example, D. M. Young, R. T. Gregory, A Survey of Numerical Mathematics (Dover, New York, 1988), Vol. I, Sec. 7.3.

For a discussion of the use of higher derivatives in numerical computations see R. W. Hamming, Numerical Methods for Scientists and Engineers (Dover, New York, 1986), Chap. 17. Also note that one common way of specifying an index distribution is in the form of a Taylor series. In the most general case, this series can be written as n2(x, y, z)=∑i=0mi∑j=0mj∑k=0mkNijkxiyjzk. In this form, it is less costly to determine higher derivatives of the index than the index itself.

For example, with 18 extra function evaluations required for finite differencing, a simple index function whose evaluation requires five or six operations (if the index is represented in the form of the series given in Note 19, only a couple of terms require this many operations) would make the number of operations for implementing the specific form of finite differencing discussed here and the differential Runge–Kutta scheme approximately equal. With 10 operations required for evaluating the index (or one of its derivatives), the differential Runge–Kutta scheme would require fewer operations.

For the mechanics of differential ray tracing through homogeneous media see, for example, A. Cox, A System of Optical Design (Focal Press, London, 1964), pp. 112–121.

E. W. Marchand, Gradient Index Optics (Academic, New York, 1978), Sec. 6.2.

A general discussion of differential ray tracing can be found in H. A. Buchdahl, An Introduction to Hamiltonian Optics (Dover, New York, 1993), Sec. 10.

See, for example, R. Kingslake, Lens Design Fundamentals (Academic, Orlando, Fla., 1978), Chap. 10, Sec. I.

Other coordinates may be chosen, and the principles presented in Subsection 2.D can still be applied. The use of Cartesian coordinates and optical direction cosines to specify position and direction is convenient for the particular method for numerical integration used here.

More generally, the matrix ∂v/∂w, where v= (v1, v2, ⋯ , vj) and w=(w1, w2, ⋯ , wk), is shorthand for the j×k matrix given by ∂v∂w=∂v1∂w1∂v1∂w2⋯∂v1∂wk∂v2∂w1∂v2∂w2⋮⋮⋱∂vj∂w1⋯∂vj∂wk.

The outer product of two vectors, v=(v1, v2, ⋯ , vj) and w=(w1, w2, ⋯ , wk), is the j×k matrix given by v⊗w=v1w1v1w2⋯v1wkv2w1v2w2⋮⋱⋮vjw1⋯vjwk.

Specifically, J is given by J=[nα′ cos(Iα′)-nα cos(Iα)]1+∂xα∂yα21/2, where Iα is the angle of incidence and Iα′ is the angle of refraction.

See, for example, J. L. Synge, Geometrical Optics—An Introduction to Hamilton's Method, Vol. XXXVII of Cambridge Tracts in Mathematics and Mathematical Physics, F. Smithies, J. A. Todd, eds. (Cambridge U. Press, London, 1962), Chap. 5, Sec. 20.

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (6)

Fig. 1
Fig. 1

Schematic representation of refraction at surface α. Quantities specified after refraction are distinguished by primes. The refractive indices before and after the surface are denoted nα and nα, respectively, and neither is assumed to be homogeneous. The point of intersection of the ray with the surface is written as (xα, yα). The optical direction cosines at this point are (pα, qα) before refraction and (pα, qα) after refraction.

Fig. 2
Fig. 2

Schematic representation of the base ray within the inhomogeneous medium. At each of the internal Runge–Kutta points a plane that is locally normal to the base ray is used for reference.

Fig. 3
Fig. 3

Relative errors in the computed values of (yβ/yα), (yβ/qyα), (qyβ/yα), and (qyβ/qyα) plotted as a function of the step parameter Δt. Since this particular Runge–Kutta scheme is globally fourth order, the error falls by 4 orders of magnitude when Δt is reduced by 1 order of magnitude. Note that the final base plane is near the focal plane [so the value of (yβ/yα) is small]; this accounts for the larger relative error in (yβ/yα).

Fig. 4
Fig. 4

Three physically interpretable quantities that completely characterize the first-order properties (with respect to the base ray coinciding with the axis of the system): the locations of (i) the front point, (ii) the rear focal point, and (iii) the rear focal length.

Fig. 5
Fig. 5

Seven physically interpretable quantities [labeled (i)–(vii)] that can be used to characterize first-order properties. These quantities are described in the text.

Fig. 6
Fig. 6

Three physically interpretable quantities [labeled (viii)–(x)]. They can be combined with the seven quantities illustrated in Fig. 5 to form a set of 10 physically interpretable quantities that completely characterize the first-order properties of this system with respect to the skew base ray used in this example.

Equations (102)

Equations on this page are rendered with MathJax. Learn more.

p2=n2[x(y), y]-q2,
M(α,β)=yβyαyβqαqβyαqβqα,
yβyα=yβyαyβzαzβzαzβzα.
xβyα=xβyβ yβyα.
pβyα=1pβ 12 nβ2xβ xβyβ+nβ2yβ yβyα-qβ qβyα,
M(α,β)=Rβ(Tβ-1Rβ-1)(Tα+2Rα+2)(Tα+1Rα+1)Tα.
Tα=yα+1yαyα+1qαqαyαqαqα.
qα/yα=0,qα/qα=1,
yα+1=yα+Lqα,xα+1=xα+Lpα,
yα+1yα=1+qαLyα,
xα+1yα=xαyα+pα Lyα,
Nα+1=1, -xα+1yα+1.
xα+1yα, yα+1yαNα+1=0,
xα+1zα, yα+1zαNα+1=0.
Lyα=1(pα, qα)Nα+1 xα+1yα+1-xαyα.
yα+1yα=1+1(pα, qα)Nα+1 qαxα+1yα+1-xαyα.
yα+1qα=qαdLdqα+L1,
xα+1qα=pα Lqα+L pαqα.
dLdqα=L(pα, qα)Nα+1 xα+1yα+1+qαpα.
yα+1qα=L1+L(pα, qα)Nα+1 qαxα+1yα+1+qαpα.
Rα=10qαyαqαqα.
qα=qα-Jxαyα,pα=pα+J,
qαyα=-xαyαJyα-J 2xαyα2,
pαyα=pαyα+Jyα,
pαyα=12pα nα2xα xαyα+nα2yα
2xαyα2=2xαyα22xαyαzα2xαyαzα2xαzα2.
pα pαyα+qα+1 qαyα=12 nα2xα xαyα+nα2yα.
Jyα=1(pα, qα)Nα 12 nα2xα xαyα+nα2yα-nα2xα xαyα+nα2yα pαpα+Jqα 2xαyα2,
qαqα=1-xαyαJqα,
Jqα=1(pα, qα)Nα qα pαpα-qα.
d2Rdt2=12 n2(R)R.
R=xy,T=dx/dtdy/dt=pq,
D(R)=12 n2(R)x,n2(R)y.
R(k+1)=R(k)+ΔtT(k)+j=1mujAj,
T(k+1)=T(k)+j=1mvjAj,
A1=ΔtD(R(k)),
Aj=ΔtDR(k)+Δtw0,jT(k)+i=1j-1wi,jAi,
j=2, 3, , m.
R(k+1)=R(k)+Δt[T(k)+16(A1+2A2)],
T(k+1)=T(k)+16(A1+4A2+A3);
A1=ΔtD(R(k)),
A2=ΔtD[R(k)+Δt(12T(k)+18A1)],
A3=ΔtD[R(k)+Δt(T(k)+12A2)].
Rα+1=R(l)+ΔtτT(l)+τ22 1-τ+τ23A1+τ33 (2-τ)A2+τ36 (τ-1)A3.
Tα+1=T(l)+b1A1+b2A2+b3A3+b4A4.
A4=ΔtD(R(l)+c0ΔtT(l)+1/2c02ΔtA3),
c0=34-14 9+16ττ-31/2,
b1=112(2c0-3)(τ-3)[c0(4τ2-9τ+6)-3τ(τ-1)2],
b2=2τ3(c0-1)(3τ3-20τ2+39τ-24)3[c0(4τ-6)+(3τ2-10τ+9)],
b3=τ2[c0(4τ-3)+τ(2-3τ)]6(c0-1),
b4=112τ(τ-1)(τ-3)2(2c0-3).
Tα=S(l)S(l-1)S(l-2)S(0),
x(0)y(0)=xαyα,
x(k)y(k)=-q(k)p(k)for(1kl),
x(l+1)y(l+1)=xα+1yα+1.
y(k+1)y(k)=y(k+1)y(k)t+y(k+1)tty(k)=y(k+1)y(k)t+q(k+1)ty(k),
y(k+1)q(k)=y(k+1)q(k)t+y(k+1)ttq(k)=y(k+1)q(k)t+q(k+1)tq(k),
q(k+1)y(k)=q(k+1)y(k)t+q(k+1)tty(k)=q(k+1)y(k)t+12 n(k+1)2y(k)ty(k),
q(k+1)q(k)=q(k+1)q(k)t+q(k+1)ttq(k)=q(k+1)q(k)t+12 n(k+1)2y(k)tq(k).
x(k+1)y(k)=x(k+1)y(k)t+p(k+1) ty(k).
N(k+1)=(p(k+1),q(k+1))0kl-1(1,-xα+1/yα+1)k=l.
[x(k+1)/y(k), y(k+1)/y(k)]N(k+1)=0,
[x(k+1)/z(k), y(k+1)/z(k)]N(k+1)=0.
ty(k)=-1n(k+1)2 p(k+1) x(k+1)y(k)t+q(k+1) y(k+1)y(k)t0kl-1-1(pα+1, qα+1)N(k+1) xα+1y(k)t-xα+1yα+1 yα+1y(k)tk=l,
R(k+1)=R(k)+ΔtT(k)+j=1mujAj,
T(k+1)=T(k)+j=1m vjAj,
R(k+1)y(k)=R(k)y(k)+ΔtT(k)y(k)+j=1m uj Ajy(k),
R(k+1)q(k)=ΔtT(k)q(k)+j=1m uj Ajq(k),
T(k+1)y(k)=T(k)y(k)+j=1m vj Ajy(k),
T(k+1)q(k)=T(k)q(k)+j=1m vj Ajq(k),
R(k)y(k)=x(k)/y(k)1,
T(k)y(k)=12p(k) n(k)2x(k) x(k)y(k)+n(k)2y(k)0,
T(k)q(k)=[-q(k)/p(k)1],
dA1dy(k)=Δt D(R(k))R(k) R(k)y(k),
D(R)R=12 2n2(R)x22n2(R)xy2n2(R)xz2n2(R)xy2n2(R)y22n2(R)yz2n2(R)xz2n2(R)yz2n2(R)z2.
Rj :=R(k)+Δtw0, jT(k)+i=1j-1 wi, jAi.
Ajy(k)=D(Rj)Rj Rjy(k),Ajq(k)=D(Rj)Rj Rjq(k),
Rjy(k)=R(k)y(k)+Δtw0, j T(k)y(k)+i=1j-1 wi, j Aiy(k),
Rjq(k)=Δtw0, j T(k)q(k)+i=1j-1 wi, j Aiq(k).
Rα+1y(l)t=R(l)y(l)+Δtτ T(l)y(l)+τ22×1-τ+τ23 A1y(l)+τ33 (2-τ) A2y(l)+τ36 (τ-1) A3y(l),
A1y(l)=Δt D(R(l))R(l) -q(l)/p(l)1,
A2y(l)=Δt DR(l)+Δt12 T(l)+18 A1R×R(l)y(l)+Δt12 T(l)y(l)+18 A1y(l),
A3y(l)=Δt DR(l)+ΔtT(l)+12 A2R×R(l)y(l)+ΔtT(l)y(l)+12 A2y(l).
Rα(b)=xα(b)(yα(b))yα(b),
Tα(b)=[nα2(xα(b), yα(b))-(qα(b))2]1/2qα(b).
Δyα(1)=(Δα(1), 0),Δqα(1)=(0, 0),
Δyα(2)=(0, Δα(2)),Δqα(2)=(0, 0),
Δyα(3)=(0, 0),Δqα(3)=(Δα(3), 0),
Δyα(4)=(0, 0),Δqα(4)=(0, Δα(4)),
Rα(j)=xα(yα+Δyα(j))yα+Δyα(j),
Tα(j)=[nα2(Rα(j))-(qα+Δqα(j))2]1/2qα+Δqα(j)
Δyβ(j) :=yβ(j)-yβ(b),Δqβ(j) :=qβ(j)-qβ(b).
M(α,β)Δyβ(1)Δα(1)Δyβ(2)Δα(2)Δyβ(3)Δα(3)Δyβ(4)Δα(4)Δqβ(1)Δα(1)Δqβ(2)Δα(2)Δqβ(3)Δα(3)Δqβ(4)Δα(4)
Δyβ(1)Δα(1)=Δyβ(1)/Δα(1)Δzβ(1)/Δα(1).
M(α,β):=yβyαyβqαqβyαqβqα0.000558998461,12.500546mm1-0.079949695 mm-111.04679141,
xαyα0.062220206mm1.0mm-1.0 mm,
pαqα0.90553851-0.3-0.3.
yβ-4.0498548mm-4.1560984mm,pβqβ0.89603525-0.37758140-0.23356611.
M(α,β)0.00662011560.005356032315.205875mm1.3921343mm0.00519148450.0156018521.4052585mm15.307419mm-0.067217237mm-10.0074985286mm-10.90660313-0.0710129360.0050730740mm-1-0.063663769mm-1-0.0609770331.0023253.
vw=v1w1v1w2v1wkv2w1v2w2vjw1vjwk.
vw=v1w1v1w2v1wkv2w1v2w2vjw1vjwk.
J=[nα cos(Iα)-nα cos(Iα)]1+xαyα21/2,

Metrics