Journal of Optometry (2014) xxx, xxx-xxx

journaj Optometry

www.journalofoptometry.org

ORIGINAL ARTICLE

Third-order aberrations in GRIN crystalline lens: A new method based on axial and field rays

Arturo Díaz del Río, Carlos Gómez-Reino, M. Teresa Flores-Arias *

Grupo de Microoptica y Óptica GRIN, Departamento de Física Aplicada (área de Óptica), Facultade de Óptica e Optometría (Campus Vida), 15782 Santiago de Compostela, Spain

Received 23 December 2013 ; received in revised form 13 August 2014

KEYWORDS

Human eye; Crystalline; Third-order aberrations; Axial and field rays

Abstract This paper presents a new procedure for calculating the third-order aberration of gradient-index (GRIN) lenses that combines an iterative numerical method with the Hamiltonian theory of aberrations in terms of two paraxial rays with boundary conditions on general curved end surfaces and, as a second algebraic step has been presented. Application of this new method to a GRIN human is analyzed in the framework of the bi-elliptical model. The different third-order aberrations are determined, except those that need for their calculation skew rays, because the study is made only for meridional rays.

© 2013 Spanish General Council of Optometry. Published by Elsevier España, S.L.U. All rights reserved.

PALABRAS CLAVE

Ojo humano; Cristalino; Aberraciones de tercer orden; Rayos axiales y de campo

Aberraciones de tercer orden en los cristalinos con gradiente de índice: nuevo método basado en rayos axiales y de campo

Resumen Este documento presenta un nuevo procedimiento para el cálculo de las aberraciones de tercer orden en los cristalinos con gradiente de índice (GRIN), que combina un método numérico iterativo con la teoría de Hamilton sobre aberraciones, en términos de dos rayos parax-iales con condiciones de contorno sobre superficies generales con límite curvado, y que se ha presentado como segundo paso algebraico. Se analiza la aplicación de este nuevo método al GRIN humano en el marco de un modelo bi-elíptico. Se determinan las diferentes aberraciones de tercer orden, excepto aquellas cuyo cálculo precisa de rayos inclinados, ya que el estudio se realiza únicamente para rayos meridionales.

© 2013 Spanish General Council of Optometry. Publicado por Elsevier España, S.L.U. Todos los derechos reservados.

* Corresponding author at: Universidade de Santiago de Compostela, Applied Physics, Facultade de Optica y Optometría, Campus Vida, Santiago de Compostela 15782, Spain. Tel.: +34 881813502; fax: +34 881813642. E-mail address: maite.flores@usc.es (M.T. Flores-Arias) .

http://dx.doi.org/10.1016Zj.optom.2014.09.003

1888-4296/© 2013 Spanish General Council of Optometry. Published by Elsevier España, S.L.U. All rights reserved.

Introduction

The most common models of the human crystalline lens are those found in the standard schematic eyes and they represent the optical structure by spherical and/or aspheri-cal end surfaces and constant refractive index. Schematic eye models are used to estimate basic properties of the eye including paraxial properties, ocular aberrations, etc. The most successful eye model was proposed by Gullstrand1 and updated by Le Grand.2 This model reproduces the Gaussian properties of an average eye. All those models have the weakness that they do not use a variable refractive index and therefore do not accurately represent the optical structure of the lens. To understand the effect of inhomogeneity of the refractive index within the lens and then on the eye as a whole, we need understand the role of the gradient-index (GRIN) nature in optical quality of lens. One important problem is that the exact distribution of the refractive index of the human lens is not well known yet. Navarro et al. proposed a GRIN lens model with concentric iso-indicial contours mimicking the external conic surfaces of the lens. The GRIN spatial distribution includes the age dependence.3 Goncharov and Dainty4 used a different approach with a fourth-order polynomial for describing the GRIN lens of a wide-field schematic eye model. Diaz et al. used a combination of polynomials and trigonometric functions for describing the refractive index distribution.5 These last two models also analyze the dependence on the age. Baharami and Goncharov6 present a new class of GRIN lens based on experimental data. The model allows, in some cases, analytical paraxial ray tracing.

It is clear that the GRIN structure may play a primary role in paraxial domain and aberrations of different orders since it permits predict real behavior of the lens. Some schematic models have shown the crucial role of aspheric surfaces in keeping aberrations.7 In the GRIN modeling of crystalline lens of the human eye, two different models are used. In one, refractive index profile is represented by a finite and discrete set of concentric shells, with a constant refractive index in each shell.8"9 In the other model, the refractive index profile is described by continuous isoindi-cial surfaces of different shapes.7,8,10"12 Some authors use numerical methods for ray tracing for evaluating radial, spherical or cylindrical refractive index gradient, and some consist of polynomial expansions of the ray trajectory inside the medium.12 The most schematic GRIN lens models had aspherical-like refracting end surfaces as expected in real lenses.13

On the other hand, the aberrations provided by an optical system in rays traversing through it may be evaluated by either a ray tracing technique or by analytical method.14 The first one involves step-by-step integration of the ray equation for rays starting from the object point and finding the intersection point in the image plane.3,12,14,15 The second one determines algebraically the aberrations of the system by equations which describe the propagation of rays through it.14,16,17 Regarding the commercial software, to our knowledge some ray tracing software (as Zemax or Code-V) does not correctly calculate the paraxial ray tracing when the optical system dealt with GRIN media. Hence, the slope and the height of rays are wrongly traced to the posterior surface of the lens. In our method we use an iterative

calculation to exactly determine the position (height) and slope of the rays at the back surface of the lens. The analytical method for determining the third-order monochromatic aberration of a GRIN system has been developed by Luneb-urg and Buchdhal.18"21 Both assumed that the system was symmetric and that the refractive index varies continuously throughout the system, that is, no interfaces between distinct media are allowed. Basing his analysis on the theory of quasi-invariants developed by Buchdahl, Sands16 has made a study of the third-order aberration for any symmetric system with any number of interfaces between a pair of inho-mogeneous media or homogeneous/inhomogeneous media, that is, for a system where no plane surfaces between media are allowed. Moore,15,22 based on the third-order aberration coefficients developed by Sands, used glasses with continuously varying refractive index in the design and construction of GRIN singlets for imaging systems.23 Thya-garajan and Ghatak17 have extended Luneburg's results on Hamiltonian theory of aberrations to the third-order aberration of inhomogeneous lenses. Both gave explicit expressions for third-order aberration in terms of two paraxial rays obeying certain specific boundary conditions on the object plane and any other reference plane of the system. These results have been applied to a GRIN rod, GRIN medium with cylindrical symmetry, limited by plane parallel end surfaces.14 Third-order aberrations in GRIN lenses with curved end surfaces may be studied in terms of two paraxial rays with boundary conditions on these non-planar surfaces.

Bahrami and Goncharov6 use a description for aberration coefficients of a thin homogeneous layer for a general GRIN lens description. Díaz et al. use the ABCD matrix for obtaining the ray tracing that allows to calculate the third-order aberrations.12,24 Recently, in 2012, Díaz et al.25 reported additional insights considering the changes that the eye suffers with age. In this paper, they evaluated separately the role of the cornea as well as the lens. For the crystalline, third-order aberrations due to the refraction in the end surfaces and to its GRIN nature have been determined. Navarro et al.7 developed an analyzed a method to obtain schematic models of individual eyes that were able to reproduce their monochromatic wave aberrations. In 2007,26 these researchers proposed and analyzed an aging and accommodating human lens plausible in terms of shape, GRIN structure and optical performance. They apply this GRIN model to develop a schematic but realistic optical model of the human lens.

Chromatic aberrations have also been analyzed in the literature.27"28 Goncharov and Dainty4 presented a mathematical method to construct a GRIN lens with iso-indicial contour following the optical surfaces with a given aspheric-ity. In this paper, the role of the GRIN structure in relation to the lens paradox is analyzed. Ocular wavefront aberrations as well as the effect of the aging on the anatomical structure of the age were studied. In 2012, Bahrami and Goncharov29 analyzed the dispersion throughout the GRIN structure of the crystalline. They use the geometry-invariant GRIN lens monochromatic model and introduce wavelength dependence of the refractive index. They developed a paraxial ray tracing method.

The aim of this paper, as a first step, is to present a novel procedure for calculating the third-order aberration of GRIN lenses that combines an iterative numerical method with the

Third-order aberrations in gradient-index crystalline lenses

Hamiltonian theory of aberrations in terms of two paraxial rays with boundary conditions on curved end surfaces and, as a second step, to apply this procedure for evaluating the algebraic form of third-order aberration coefficients of the human lens in the framework of that asymmetric bi-elliptical continuous GRIN model that provides a close simulation to the real one. Let us underline that the method proposed in this work is completely general in terms of the back and front surfaces shape of the human lens, because no restrictions have been imposed on them. Moreover, we can determine the third order aberration for any input point on the front surface and any entrance angle. In this proposal, both, refraction and propagation of light throughout a GRIN lens described by a continuous refractive index variation are considered. In the particular case analyzed in section 'A particular case application: the GRIN human lens', back and front elliptical end surfaces for the human lens have been assumed.

Determination of the third order aberrations

If we assume that the lens is rotationally symmetric about the z axis, the refractive index distribution can be expressed

Ho = H(0, 0); Hi = —

that is

The Hamiltons's equations are

dx „ 9H

— = 2p— dz r dv

dp „ 9H

-Ç- = —2x— dz 9u

It is not possible to solve the equations exactly and it will be necessary to derive approximate solutions. The lowest order approximation is the paraxial approximation and the deviation from this approximation determines different orders of aberrations.

To calculate the third-order aberration we expand H in terms of u and v up to quadratic values to get

H = H =+HiU + H2V + ^(tfiiu2 + 2Hnuv + H22V2) + • • • (7) where H,'s and H,/s are given by14

92H 92n

"du2" u=0,v=0 — 9Û2 u=0

92H 1

"dv2 u=0,v=0 4n0(z)

_ 9n "— dU

H2 = dH 2 dv

u=0,v=0

2no(z)

2n2(z) du

in general form as obey a relation of the form22

n(u, z) = no(z) +ni(z)u + n2(z)u2 + •• • = ni(z)U (1)

where ni(z) are polynomials in z as follows

n, (z) = ]Tny (z)zf (2)

u = x2 + y2 (3)

The expressions for ni depend on the optical modeling of inhomogeneity of the refractive index within the lens. n0(z) is the axial refractive index distribution, n1 (z) affects the refraction of paraxial rays, n2(z) affects the refraction of non-paraxial rays to the third-order aberration and remainder ni affect the higher-order aberrations.

Following Luneburg, the Hamiltonian of the lens is represented by18

To determine one ray out of all the possible rays that can pass through the lens, we need know the value of (x, y, p, q) at two surfaces, the end surfaces of the lens, as a set of boundary conditions for solving the ray equation.

Hence, we expand the position and the optical direction cosine of the ray in ascending power of coordinates of the two points on the surfaces

P = Pi + P3 + P5 +----with p = x or y

l = l1 + /3 + /5 +----with l =p or q

(9a) (9b)

H = v/n2(u,z) — v

where v=p2 + q2, p and q representing the optical direction cosines of the rays along the x and y directions respectively at the point (x, y, z),

where subscripts represent the order of the term, p1 and l1 represent linear terms (paraxial approximations) in x(z0), y(z0), x(z2) and y(z2), (see Fig. 1), p3 and l3 cubic terms (third-order aberration) in the same variables and so on. Substituting Eqs. (9a), (9b) and (7) in (6a) we get

/1 + /3 + • • • = 2(l1 + /3 + • • -)(H2 + H12P2 + H22/2) (10)

where dot represents derivative with respect to z Then,

/1 = 2H2/1 (11a)

/3 = 2[H2/3 + (H12p2 +H22/2H1] (11b)

Similarly, using Eq. (6b), we get /1 = —2H1P1 (12a)

/3 = —2[H1P3 + (Hhp2 + H12/1 )p1 ] (12b)

u 0,v 0

IP3b(Z2) Plb(Z2>

Figure 2 Scheme of variables used to evaluate the third-order aberrations.

Figure 1 Evolution of the axial and field rays in a bi-elliptical GRIN medium.

We derive expressions for the third-order aberrations in terms of two paraxial rays obeying certain specific boundary conditions. Not skew rays have been used.30 These rays are solutions of Eqs. (11a) and (12a) with the following boundary conditions (see Fig. 1 )

(i) Front part:

Axial ray hf (zo) = 0; hf (01) = 1 Field ray Hf (zo) = bo; Hf (01) = 0;

(ii) Back part

Axial ray hb(01) = 0; hb(z2) = ¿2

Field ray ^(01) = 1; Htfe) = 0;

(13a) (13b)

(14a) (14b)

Providing the lens divided into front and back parts and limited by two curved surfaces joined smoothly at the equator.

The solutions under these boundary conditions are represented by ray position h(z) and optical direction cosine &(z) for the axial ray and by H(z) and Q(z) for the field ray, respectively (see Fig. 1). Axial and field rays are two linearly independent solutions of the paraxial ray equation, and any other paraxial ray, in the two parts, is given by a linear combination of these rays, that is

P1/(z) = Hf (z) (z)| /1/(z) = 0/ (z) + M/ (z) j for the front part (0 < z < a1 ) and P1b(z) = hb(z) +b1H6(z)] /1b(z) = ^ (z) + bA(z) I

(zo,bo) ^(01, M (15)

(a1,b1) ^(z2,b2) (16)

for back part (a1 < z < a2).

To evaluate third-order aberrations we have carry out the following steps (see Fig. 2)

1. Solve the Hamilton equations for paraxial domain p1f and

l1f and for third-order aberration p3f and l3f in the front part of the lens.

2. Evaluate p1b and l1b in the back part of the lens for any ray from pf(a) = P1f(a) + P3f(a1) at the equatorial plane and lf (a1) = l1f(a1) + l3f(a1).

3. Evaluate p3b and l3b in the back part of the lens for any ray from pf(a1) at the equatorial plane, previous calculation of z2.

4. Find the ray high pb(z2) for the third-order aberration on the back surface of the lens by an iterative numerical method.

To evaluate third-order aberration in the front part of the lens, we use a procedure similar to the one used by Luneburg18 in third-order aberration theory and obtain

[P3^/ - //] = 2[(h12p2/ +H22/2/)// + (hhp2/ +H12/Î/ )p1/h/ ]

Integrating the above equation from z = z0 to z = a1 and using the fact hf (z0) = 0; hf(a1) = 1 and p3f(z0) = 0 we get

P3^f - hf/3f 101 = p3f(01 (01) - l3f (01)

= ^[(H12P?/ + H22 /2/ /

+ (H11P2/ +H12/?/ )p1/h/ ]dz

Similarly

^ [P3/0/ - //] = 2[(H12P2/ +H22/2// + (H11P1/ + H12 /1/ )P1/H/ ]

Integrating from z = z0 to z = a1 and using the fact Hf (z0)= b0; Hf(a1) = 0 and l3f (z0) = 0, we have

p3f0f -Hfl3f 101 = p3f(01)0f (01) = 2 J[(H12p2f +H22/?f)f

+ (H11P2/ +H12/Î/ )P1/H/ ]dz

Third-order aberrations in gradient-index crystalline lenses

From Eq. (15) and after straightforward calculation, Eqs. (18) and (20) can be written as

p3f (0) = [A/ + B/b1 + + f/b3] (21)

l3f (fl1) = F^TT// (a1) — E/0/ (a1) + / (01) 0/ (01)

— C/0/ (01))b1 + (// (01) — 6/0/ (01))b1 + (E/fy (01)

— A/0/ (01))b?l

where the coefficients Af, Bf, Cf and Ef are, in turn, given

A/ =2 J [Huh4 + 2H12h2^2 + H22tf4]dz

6/ = 6 J [HuH/h3 + H12// (H/fy + h/0/ ) + H220/^3]dz

C/ = 2 J {3[HnH/2h2 + H2202^2] + H12[H/2^/2 + h202 z0

+ 4H/h/0/^/ ]}dz (23c)

E/ = 2 J [HuH/h/ + H12H/0/ (H/fy + h/0/ ) + H2203^/ ]dz z0

in Eqs. (21) and (22), superindex * denotes change of variables Hf by hf and 0f by fy for a given coefficient. Note that Cf = Cf*. Both equations give the third-order aberration for the ray position and the optical direction cosine at the equatorial plane. Eqs. (21) and (22) can be rewritten, taking into account second line in Eq. (15), as

P3/ (01) =

I3/ (01) =

61 (zo) 0/ (01)

62 (zp) 0/ (01)

where e1 and e2 are polynomials in third-order of l1f(z0), given by

•1(z0) = ( /0}' '' E'

C ; — 3 E

/ fy(zo) /

0/ (zq) fy (zo)

0/ (zo) fy (zo)

l1/ (zo) fy (zo)

l1/(zo) I A ; fl* 0/(zo)

fy (zo) / (zo)

6;—2/

/ (zo)

0/ (zo) fy (zo)

62(zo) = //(01) — A/0/(01)] (/p})

// (01) — 6/0/ (01) — 3 E/fy (01)

—A/0/ (01)

0/(zo) fy (zo)

l1/(zo) fy (zo)

(01) — C/0/ (01) — (01)

—6/0/(01)) ( ) + 3//(01)

—A/0/(01))

fy (zo)

0/ (zo) fy (zo)

l1/(zo) fy (zo)

+ (01)

— E/0/ (01) — / (01) — C/0/ (01))

0/(zo) fy (zo)

+(// (01) — 6/0/ (01))

0/ (zo) fy (zo)

—/ (01)—a/0/ (01)^ /q))

Then, the position pf and the optical direction cosine lf of the ray at the equatorial plane for the third-order aberrations are given by

P/ (01) = b1 +P3/ (01) l/ (01) = l1/(01) + 13/ (01)

(26a) (26b)

We evaluate now the third-order aberrations in the back part of the lens. From the procedure used in the front part of the lens to find third-order aberrations, we obtain

(p3b#b — hb/3b) 102 = p3b(z2)#b(z2) — 62/3^2) — p3b (01)^6 (01)

= ^[(H12P2b +H22l2b)l1bfyb

+ (H11P2b + H12l1b)P16hb]dz

(P3b0b — Hbl3b) |01 = P3b(z2)0b(z2) — P3b(01 )0b(01 ) + l3b(01)

= 2J [(H12P2b +H22lib)l160b + (H11P2b 01

+ H12l2b)P16Hb]dz (28)

After calculations, Eqs. (27) and (28) can be expressed as

P3b(z2)fyb(z2) — b2l3b(z2) = Ab + 6^1 + C6b2 + Ebb3 (29) P3b(z2)0b(z2) + l3b (01) = Eb + Cbb1 + 6bb2 + Abb1 (30)

where the conditions l3b (a1 ) = p3b (a1 ) = 0 have been used and the coefficients Ab, Bb, Cb and Eb are given by

Ab =2 J [Hnhb(z) + 2H12hb(z)^2 + H22tf4(z)]dz (31a)

ßb = 6J [H11H(z)hb(z) + H12hb(z)^b(z)(Hb^b + Mb)

+H22^3(z)0fc(z)]dz

Cb = 2J {3[H11Hb(z)hb(z) + ^(z^z)]

+ H12[Hb(z)^2(z) + h2(z)0b(z)

+ 4Hb (z)hb (z)^b (z)0b (z)] }dz

Eb = 2 J [H11Hb(z)hb (z) + H12Hb(z)0b(z)(Hb(z)^b(z)

+ hb(z)0b(z)) + H2203(z)#b (z)]dz

Superindex * indicates change of variable Hb by hb and 0b by $b for a given coefficient as commented in aberrations for the front part.

From Eqs. (29) and (30), it follows that the third-order aberrations for position and optical direction cosine of the ray at z = z2 are given by

1 2 3 P3b(z2) = ) [Eb + Qb + ß*b1 + Abb3]

/3b (z2) = 1 i ?4§4[Eb + Cbb1 + ßbb1 +Abb1]

b2 V 0b(z2) -[Ab + ßbb +Cbb2 +Ebb3]} Eqs. (32) and (33) also can be expressed as «1 (01)

P3b(z2) =

/3b(z2) =

0b (z2) (01 )

0b (z2) where

«1 (01 ) = Eb+Cbb1 +ßb b1 +Abb1

1 2 3 «2(01) = r {^b (z2)[Eb +Cb b1 + ßbb? + Abb?] b2

0b(z2)[Ab + ßbb1 + Cbb1 + Ebb3]}

The position pb and the optical direction cosine lb of the ray at the curved output surface z= z2 for the third-order aberrations are given by

Pb(z2) = P1 b (z2) + P3b(z2) /b (z2) = /1b(z2) + /3b(z2)

Front part

Back part

Figure 3 Structure of the refractive index of the human lens in the sagittal section considered as a bi-elliptical iso-indicial curves.

For calculating the angle of the ray out of the lens, we use the Snell law, given by

>,out = sin

sin /3

with ne the value of the refractive index at the edge of the lens and nout the refractive index of the medium in contact with this surface.

The focal plane of the lens is determined by an iterative process for a set of parallel rays in the paraxial region, as the intersection of these rays with the axial axis.

A new algorithm was developed to locate the value of the abscissa z2 (see Fig. 1), corresponding to the intersection of the third-order aberrated ray trajectory with the posterior surface of the lens.

In order to determine the third-order aberrations we have used the definitions and formula given by Marchand.31

A particular case application: the GRIN human lens

In order to evaluate the third-order aberrations in a particular case, we analyze in this section a human lens in the framework of the asymmetric bi-elliptical continuous GRIN as proposed by Gomez-Reino et al.32 The aperture stop was supposed to be in contact with the lens. The lens is considered as a whole, delimited by two elliptical surfaces, frontal and back. Because of the partition of the lens, we have introduced an imaginary plane surface between the anterior and posterior parts to demarcate the start of the posterior section (see Fig. 3). This surface has infinite radius of curvature and no asphericity. It has not been considered the dependence of the crystalline on the age and we have studied a case accommodation free.

Third-order aberrations in gradient-index crystalline lenses 7

Figure 4 (a) Longitudinal spherical aberration (LSA) and (b) transversal spherical aberration (TSA) versus the input angle Calculations have been made for the coefficients c0 = 1.406, c1 = -0.0201416, c2 =0.0001423, c3 = -0.0000007 and c4=0, ai =1.7mm and a2 = 1.9 mm and a lens semiaperture of 4.5 mm.

Assuming that the human lens has symmetry of revolution, the Eq. (2) can be rewritten as

n(u, z) = n0nr where

n0 = cjZJ = C0 + C1Z + C1Z2 + •

being cj the coefficients of the power series and

nr = 1 - 2+ 2 ^4u2

with u = r2, g the gradient parameter which characterizes the refractive index distribution and p the coefficient representing the shape of the refractive index distribution.11

Taking Eqs. (41)-(43) into account, with the help of a developed Matlab iterative program, we can determine the trajectory of any meridional ray that impinges at the front part of the crystalline using a power series. The convergence of the series is important, because the aberrations are calculated from the difference of two ray heights with accuracy higher than 10-6.33

Once the path of a ray affected by aberration is calculated, we can determine the value of the third-order aberrations using the expressions determined by Marchand.31 At this point is important to notice that in this work, we do not use skew rays and therefore, the sagittal curvature and the astigmatism could not be approach, this is also the reason because there is not D coefficient in Eqs. (21)-(22) and (32)-(33).

For performing the aberrations for the particular case of the GRIN human lens, values of central refractive index of 1.406 and edge index of 1.386 were taken from the Gull-strand schematic eye. The coefficients cj of the power series in the refractive index profiles are taken to be c0 = 1.406, c1 = -0.0201416, c2 =0.0001423, c3 = -0.0000007 and c4 =0. The radius of the crystalline is taken to be 4.5 mm and the semiaxes along the z-axis of the asymmetric bi-elliptical iso-indicial curves a1 and a2 have the values of 1.7 and 1.9 mm,

respectively. The values of the gradient parameter are taken to be g = 0.03775 mm-1 and p = 3/4.32

Fig. 4a shows how the longitudinal spherical aberration (LSA) increases, in absolute value, with the input angle § (defined as the angle between the impinging ray and the axial axis, see Fig. 2). For the paraxial region, values up to -0.1 mm have been obtained. For angles outside of the paraxial region, the dependence of the LSA with the input angle is nearly parabolic. Note that values for the LSA are negatives due to the GRIN nature of the crystalline, as reported in Ref.34 Fig. 4b depicts a similar behavior for the transversal spherical aberration (TSA), however, in this case, the TSA is nearly constant in the paraxial region, and the value of the TSA for 0.1 angle is -18 x 10-5 mm. For angles between 0.1° and 0.4°, its increase is quasi-linear and beyond this value a higher increase with § is observed, showing a parabolic behavior. The coma has been calculated at the fovea, assuming this to be 5° from the optical axis, and with pupil centered on the optical axis. Fig. 5 shows the coma versus p0. The coma increases quasi-linearly as the height of the input ray increases, it means as we move away from the paraxial area, achieving values up to 0.021 mm for a 2.1 mm input ray height.

In order to compare these values with those of other author, we can express the ocular aberrations by means of Zernike coefficients using the relationship's25

C0_ SI °4 = 6V5

°3 = 6V2

(44a) (44b)

being C4 and C] the fourth-order spherical and lateral (horizontal coma), respectively and SI and SII the Seidel sums for spherical aberration and coma, respectively.

For determining SI and SII from the values presented in Figs. 4a, b and 5, we use expressions and definitions tacked from pages 502-504 of the book Optical Shop Testing.35

In the case of the spherical aberrations, for example for a value of an input angle of § = 0.66°, LSA is equal to -8.8 mm and TSA is equal to -0.1 mm. In this case a value

Figure 5 Comma versus the height of the input ray a,. Cal- Figure 7 Meridional curvature versus the input angle tf. Calculations have been made for the same parameter that Fig. 4, culations have been made for parameters of Fig. 4. for the fovea, assuming it is at 5°.

of C, = -0.169 ^m is achieved. For the interval between tf = 0.35°-0.80°, with the values of Fig. 4, C, takes values in the interval -0.0031 to -0.499 ^m. These values are in good agreement with those presented by other authors.4,25 In the case of the coma, for the interval between 0.8 and 2 mm of the input height (p0), the values obtained expressed in the standard representation (C]) are in the interval -0.0047 to -0.0536 ^m. These values are lower than those given by other authors.4,25

The distortion is presented in Fig. 6. It has been calculated from values of tf between 0° and 20°. Dependence quasi-parabolic has been observed. Distortion increases as we separate from the paraxial region.

The meridional curvature (Zm) has been calculated according the Eq. (7.3) given at the reference.31 Fig. 7shows a parabolic behavior of the meridional curvature with the input angle tf. Values of the Zm vary in the range -5 to 5 mm for input angles in the interval 1°-21°. The sagittal curvature

Figure 6 Distortion versus the input angle tf. Calculations have been made for parameters of Fig. 4.

(Zs) requires the consideration of a skew ray, so we cannot calculate it, as explained in a previous paragraph.

To our knowledge, this is the first time that a general analytical method for determining the human lens aberrations, from its GRIN nature in the framework of the asymmetric bi-elliptical iso-indicial model, without restriction on the back and front shape surfaces has been presented.

Conclusions

A novel procedure for calculating the third-order aberration of GRIN lenses that combines an iterative numerical method with the Hamiltonian theory of aberrations in terms of two paraxial rays with boundary conditions on a general curved end surfaces and, as a second algebraic step has been presented. Application of this new method to a GRIN human lens has been made in the framework of the bi-elliptical model. The different third-order aberrations have been determined, except those that need for their calculation skew rays, because the study has been made only for meridional rays. The particular case of the crystalline human lens has been analyzed. The method proposed presents the advantage of considering a continuous GRIN media with any kind of end surface (for our case applied to aspheric curved ones). We can determine the output position and slope of any ray that travel throughout the human lens. Refraction on back and front end surfaces is considered. Difficulties of ray tracing with GRIN continuous lens have been overcome thanks to this combined analytical and iterative method that allows exactly determine the height and slope of any ray traveling throughout the lens.

The LSA presents parabolic increases as the input angle increases, as well as the TSA does. The coma increases quasi linearly with the height of the input ray at the front of the crystalline. The distortion and the meridional curvature increase smoothly with the angle at the front of the human lens.

Third-order aberrations in gradient-index crystalline lenses

Conflicts of interest

The authors have no conflicts of interest to declare.

Acknowledgment

This work has been partially supported by the project

EM2012/019 Xunta/Feder Spain.

References

1. Gullstrand A. The optical system of the eye. In: von Helmholtz H, ed. Handbuch der Physiologische Optik. vol. 1, 3rd ed. Hamburg: Voss; 1909:350-358.

2. Le Grand Y,El Hage SG. Physiological Optics. New York: Springer-Verlag; 1980.

3. Navarro R, Palos F, González L. Adaptive model of the gradient index of the human lens. I. Formulation and model of aging ex vivo lenses. J Opt Soc Am A. 2007;24:2175-2185.

4. Goncharov AV, Dainty C. Wide-field schematic eye models with gradient-index lens. J Opt Soc Am A. 2007;24:2157-2174.

5. Díaz Ja, Pizarro C, Arasa P. Single dispersive gradient-index profile for the aging human lens. J Opt Soc Am A. 2008;25:250-261.

6. Bahrami M, Goncharov V. Geometry-invariant gradient refractive index lens: analytical ray tracing. J Biomed Opt. 2012;17, 055001-1(1-9).

7. Navarro R, González L, Hernández-Matamoros JL. On the prediction for optical aberrations by personalized eye models. Optom Vis Sci. 2006;83:371-381.

8. Smith G. The optical properties of the crystalline lens and their significance. Clin Exp Optom. 2003;86:3-18, references therein.

9. Masajada AP. Numerical study of the influence of the shell structure of the crystalline lens on the refractive properties of the human eye. Ophthal Physiol Opt. 1999;19:41-49.

10. Pérez MV, Bao C, Flores-Arias MT, Rama MA, Gómez-Reino C. Gradient parameter and axial and field rays in the gradientindex crystalline lens models. J Opt A: Pure Appl Opt. 2003;5:S 293 -S 297.

11. Blaker JW. Toward an adaptive model of the human eye. J Opt Soc Am. 1980;70:220-223.

12. Díaz JA. ABCD matrix of the human lens gradient-index profile: applicability of the calculation methods. Appl Opt. 2008;47:195-205.

13. Smith G, Pierscionek BK, Atchison DA. The optical modelling of the human lens. Ophthal Physiol Opt. 1991;11:359-369.

14. Sodha MS, Ghatak AK. Ray tracing and aberrations in lens-like media. In: Inhomogeneous Optical Waveguides. New York and London: Plenum Press; 1977 [chapter 9].

15. Moore DT. Ray tracing in gradient-index media. J Opt Soc Am. 1975;65:451-455.

16. Sands PJ. Third-oder aberrations of inhomogeneous lenses. J Opt Soc Am. 1970;60:1436-1443.

17. Thyagarajan K, Ghatak AK. Hamiltonian theory of third-order aberrations of inhomogeneous lenses. Optik. 1976;44:329-342.

18. Luneberg RK. Mathematical Theory of Optics. Providence, Rhode Island: Class Notes for Brown University; 1944. Summer.

19. Buchdahl HA. Optical Aberration Coefficients. New York: Dover; 1969.

20. Buchdahl HA. An Introduction to Hamiltonian Optics. London/New York: Cambridge University Press; 1970.

21. Welfor Frs WT. Aberration of Optical Systems. Bristol/Philadelphia: Adam Hilger; 1989.

22. Ducan T. Moore desing of singlets with continuously varying indices of refraction. J Opt Soc Am. 1971;61:886-894.

23. Moore DT, Sands PJ. Third-order aberrations of inhomogeneous lenses with cylindrical index distributions. J Opt Soc Am. 1971;61:1195-1201.

24. Sorroche F, Díaz JA, Fernández-Dorado J, Arasa J. ABCD Matrix for calculating third-order aberrations of gradient index optical element. Proc SPIE. 2010;7652, 76522X-1(1-9).

25. Díaz JA, Fernández-Dorado J, Sorroche F. Role of the human lens gradient-index profile in the compensation of third-order ocular aberrations. J Biomed Opt. 2012;17, 075003-1(1-11).

26. Navarro R, Palos F, González LM. Adaptive model of the gradient index of the human lens. II. Optics of the accommodating aging lens. J Opt Soc Am A. 2007;24:2911 -2920.

27. Liou HL, Brennan NA. Anatomically accurate finite model eye for optical modeling. J Opt Soc Am A. 1997;14:1684-1695.

28. Díaz JA, Pizarro C, Arasa J. Single dispersive gradientindex profile for the aging human lens. J Opt Soc Am A. 2008;25:250-261.

29. Bahrami M, Goncharov AV. Geometry-invariant gradient refractive index lens: analytical ray tracing. J Biomed Opt. 2012;17:055001.

30. Marchand EW. Third-order aberrations of the photographic Wood lens. J Opt Soc Am. 1976;66:1326-1330.

31. Marchand EW. Gradient Index Optics. United Kingdom: Academic Press; 1978 [chapter 7].

32. Gomez-Reino C, Perez MV, Bao C. Gradient-Index Optics: Fundamentals and Applications. Berlin: Springer Verlag; 2002 [chapter 8].

33. Smith G. Paraxial ray tracing in gradient-index (GRIN) media: problems of convergence. J Opt Soc Am A. 1992;9:331-333.

34. de Castro A, Ortiz S, Gambra E, Siedlecki D, Marcos S. Three-dimensional reconstruction of the crystalline lens gradient index distribution from OCT imaging. Opt Express. 2010;18:21905-21947.

35. Malacara D. Optical Shop Testing. 3rd ed. New York: Wiley; 2007.