Cataract surgery is one of the most common surgical procedures nowadays, during which the patient’s crystalline lens, which has become opaque due to increased light scattering, is replaced with an artificial IntraOcular Lens (IOL). As cataracts spread across people in their lower, more active, age brackets as well, there is a growing demand for premium lenses to improve achievable image quality and to resolve focusing without eyeglasses. Diffractive IOLs offer a possible solution by simultaneously creating multiple focal points to provide sharp vision at both near and far distances. In this article, we demonstrate how to extend the capabilities of Zemax OpticStudio by using User-Defined Surface (UDS) DLLs to provide a realistic model of relief-type diffractive lenses. Finally, we also discuss the advantages of the applied zone-decomposition model over the order-decomposition method using the built-in diffractive surface types.

Authored By Csilla Timar-Fulep

Downloads

Article Attachments

The Cpp source code of the User-Defined Surface DLL can be downloaded from the Code Exchange:
DLL (User-Defined Surface): Relief-Type Diffractive Surface

Introduction

According to the increasing demand from cataract patients, IOL manufacturers invest more time and resources into the research and precise design of premium type lenses. In order to provide good image quality at a wide range of object distances, the ultimate goal would be to reproduce the accommodation capabilities of the original crystalline lens. Directly imitating the natural processes of the human eye with artificial elements faces several challenges, thus it is still an unsolved problem. However, diffractive IOLs can provide a readily available solution for multiple viewing distances at the same time. This article shows how to implement a realistic model of relief-type diffractive lenses based on the real surface shape using ray tracing together with diffraction analysis, and it demonstrates the advantages of the model for a comprehensive evaluation of the system performance.

Basic design concept

Order-decomposition

The built-in diffractive surface models in the sequential mode of OpticStudio rely on order-decomposition, where a single diffraction order needs to be selected, and then the diffractive power is represented by an additional phase contribution, independent of the refractive index and the surface sag. Using this method, the propagation of orders can be modelled either by rays from the object to the image, or by scalar diffraction from the exit pupil. This approach provides a simple solution for analyzing individual orders, which is beneficial in applications using a single target diffraction order. The working principle, and an application example of designing diffractive intraocular lenses using this representation are discussed in detail in these knowledgebase articles:
How diffractive surfaces are modeled in OpticStudio
Using diffractive surfaces to model intraocular lenses

However, there are some deficiencies of the above-mentioned order-decomposition model. First, as the phase function simply applies an extra bending to the rays on top of the base refractive or reflective surface, this model does not consider the real ray path through the diffractive element, and therefore wavelength dispersion as well as certain other aberrations are neglected. Besides, this surface model does not account for the diffraction efficiency. Finally, a multi-configuration system has to be created to model different diffraction orders one by one.

Zone-decomposition

On the contrary, using zone-decomposition, diffraction into multiple orders can be accurately considered at once, and this method inherently accounts for wavelength dispersion and diffraction efficiency by modelling the actual shape of the diffractive element. This enables the creation of advanced IOL models, where the different orders are designed to provide sharp vision for multiple viewing distances, thereby substituting accommodation.

The zone-decomposition model assumes that the zone widths of the diffractive element are much larger than the wavelength, and that the optical properties behave smoothly within the zones. In this case, geometrical optics approximation and therefore ray tracing can be used to describe the propagation from one side of the diffractive surface to the other. This also means that the zones can be regarded as conventional refractive/reflective optical elements in the near-field, while in the far-field light distribution can only be calculated by scalar diffraction analysis.

In OpticStudio, PSF calculations realize this exact process, where ray tracing is performed through the system and additionally a one-step diffraction analysis is considered from the exit pupil to the image plane. Since the phase changes through the diffractive element are calculated based on geometrical optics, zone-decomposition has the most relevance when the diffractive surface is located at the exit pupil, or at one of its conjugate positions (entrance pupil or aperture stop).

IOL design and simulation is a suitable use case, which meets the above criteria, as the implanted artificial lens is usually placed right after the pupil, which is the aperture stop of the eye. According to the common practice, when Fresnel diffraction between the aperture stop and the exit pupil is neglected, zone-decomposition can be used to efficiently model diffractive IOLs.

Diffractive surface model using UDS DLL

In order to take advantage of the above-described zone-decomposition method, we implemented a new User-Defined Surface DLL, where the sag profile of relief-type diffractive surfaces can be described analytically. Besides the accurate performance analysis of diffractive optical elements (DOEs), the parametric shape representation using the UDS DLL also enables optimization and tolerance analysis of these diffractive surfaces.  Further details about how to extend the capabilities of OpticStudio using custom DLLs and how to compile new solutions can be found in these articles:
Custom DLLs in OpticStudio: An overview of user-defined surfaces, objects, and other DLL types
How to compile a User-Defined DLL

When using sequential surface DLLs, there are 10 different ways in which OpticStudio interacts and exchanges data with the DLL. These scenarios represent general information, parameter name, and safe data transfer, as well as layout plot, paraxial and real ray trace calculations. The different functionalities are defined under different cases of the DLL.

In this model, we applied a simple rotationally symmetric diffractive structure with a uniform relief step height, added on top of a Standard surface representing the substrate. To enable model comparison with built-in OpticStudio solutions, we described the relief shape by Even Asphere polynomials. Accordingly, the surface sag is given by the following equation:

$$z=z_{subs}+z_{DOE}=\dfrac{cr^2}{1+\sqrt{1-(1+k)c^2r^2}}+mod(a_1r^2+a_2r^4+a_3r^6+a_4r^8+a_5r^{10},h)$$

In the above equation, mod represents the modulo function, c is the curvature, i.e. the reciprocal of the radius, k is the conic constant, r is the radial coordinate, and h is the uniform relief step height.

The a_i coefficients of the Even Asphere polynomials, the h step height, and the propagation algorithm indicator are first defined under Case 1, Parameter column header names, in the DLL. Then, Case 3 is used to describe surface sag based on the above equation for drawing purposes in the layout plots. Case 4 should account for paraxial ray trace results, but since diffraction analysis is required on top of ray tracing for the zone decomposition approach, which is only available for real ray trace, we ignored this step. This means that in paraxial approximation our model behaves as a Standard surface. Finally, under Case 5, the real ray trace results are computed. For this, we implemented two solutions, an approximative analytic and an iterative algorithm, which are discussed below.

Ray propagation algorithm

In case of complicated surface shapes, the ray–surface intersection coordinates cannot be determined analytically, and therefore for the built-in surface types other than the Standard surface OpticStudio applies an iterative algorithm to find a numerical solution. This can be one approach for the User-Defined DLLs as well. However, since the iterative method is computationally less efficient than the direct calculation, besides the usually applied iterative solution, we also implemented an approximative closed-form solution based on local linearization [1, 5].

In the latter alternative algorithm, we handle the sag of the substrate and the additional relief separately. First, we determine the exact ray intersection coordinates with the substrate (x₀, y₀, z₀), which can be done analytically, since the substrate has a Standard surface shape. Then, as a next step, we estimate the ray-relief intersection point based on the local relief height (Δz=z_DOE(x₀,y₀)), and the slope at the given location (x₀, y₀, z₀+ Δz). The estimated intersection point (x, y, z) with the tangent plane can be calculated analytically again by solving linear equations. This direct approximative calculation can be 30% faster than the default iterative method, without causing any significant errors in the results. The process is visualized in the image below.

Simulating the IOL performance in eye model

Bifocal IOL design

To demonstrate the applicability and advantages of the new diffractive surface DLL, we implemented an ideal diffractive lens model based on the literature [1], which equally distributes light into the zeroth and first diffractive orders at the design wavelength. (Although these two orders have equal power, the applied diffractive surface sends a small fraction of power to higher orders too.)  According to the ISO Standard for intraocular lenses, we set the primary wavelength to the e line, i.e. λ₀=546.07 nm. The IOL was designed to have a base refractive power of P₀=22.5 D with a diffractive addition of Padd=3.5 D. We modelled the lens material, Benz25, using Model material solve with refractive index of n=1.4625, while the surrounding media, natural saline, was described by Model material solve with refractive index of n₀=1.3343.

To achieve an ideal lens, where the first diffraction order is focused at a distance of EFL, the optical path difference at the boundary of zones j and j+1, relative to the spherical Gaussian reference wavefront, has to be jλ, which can be geometrically expressed as the equation below:

$$\sqrt{r_j^2+EFL^2}-EFL=j\lambda$$

This means that the zone boundaries are at the following locations:

$$r_j^2=2\cdot EFL\cdot j\lambda+(j\lambda)^2$$

Since EFL>>λ, we can neglect the latter part and apply the approximation below:

$$r_j^2 \approx 2\cdot EFL \cdot j\lambda$$

Therefore, the optical path addition of the diffractive surface to realize the above r_j zone edges and to produce 100% diffraction efficiency for the first order is as follows:

$$\Delta OPL=n\cdot mod(\dfrac{r^2}{2EFL},\lambda)$$

To realize equal diffraction efficiency for the zeroth and first orders, a constant α=0.5 multiplication factor must be taken into account as well. So, the surface sag of the ideal bifocal diffractive IOL can be described as:

$$z_{DOE}(r)=mod(\dfrac{\alpha n \cdot r^2}{2EFL(n-n_0)}, \dfrac{\alpha \lambda_0}{n-n_0})$$

Based on this, in the Lens Data Editor, we set the 2^nd order polynomial coefficient to a₁=αn/2EFL(n-n₀)=6.82E-3, and the step height to h=αλ₀/(n-n₀)=2.13E-3 mm. The sag addition of the diffractive profile, i.e. the Surface Sag map with the Base Radius being removed, is shown in the image below:

Finally, the base radius of the lens substrate was calculated by using the Lensmaker equation assuming a symmetric double-convex lens with 1.0 mm thickness and P₀=22.5 D base power. This results in a radius of 11.353 mm. The conic constant was set to k=0 on the front surface of the lens, and it was optimized on the rear surface to achieve diffraction-limited performance for the zeroth diffraction order, which resulted in the value of k=-5.8.

ISO standard eye model for production-line testing

In order to validate the IOL model, we incorporated the diffractive UDS DLL in the ISO 11979-2 standard eye model, which is designed for production-line testing of the optical properties of ophthalmic implants [2]. The model eye contains a virtually aberration-free cornea, after which the IOL is placed in a liquid medium contained between two flat windows:

The back focal length of the system was optimized for the smallest RMS Wavefront error using the Quick Focus tool, which results in an image location halfway through between the zeroth and first order focal points. This also justifies that on top of ray tracing, scalar diffraction analysis is required to properly describe the far-field behavior of the lens. In addition, at this intermediate position, the Wavefront Map clearly shows the half wave difference between the adjacent zones, which is in good agreement with the theoretical expectations.

When the diffraction effects are considered, for example via FFT PSF or MTF calculations, then all diffraction orders are accurately modelled within a single configuration, and the diffraction efficiency is inherently taken into account by the nature of the model. Results of the FFT Through Focus MTF analysis at 50 lp/mm frequency are shown below. The 0.34 peak diffraction efficiency at the zeroth and first order focal planes are close to the theoretical values.

Comparison with built-in diffraction surface models

To compare the new realistic UDS DLL model with the built-in diffractive surface models in OpticStudio, we created a multi-configuration system, where we used the Binary 2 diffractive surface type to describe the phase addition in the different orders. The phase profile corresponding to the same ideal bifocal lens as discussed before, can be described by the following equation:

$$\Phi(r)=mod(\dfrac{\pi P_{add} \cdot r^2}{\lambda_0}, 2\pi)$$

For further references, the following knowledgebase article discusses the conversion between the sag and phase profiles of diffractive elements in detail, and also provides a ZPL macro for the calculations:
How to calculate the sag of a diffractive optical element with a macro

The model contains 5 configurations corresponding to the 0^th and 1^st diffraction orders using the Binary 2 phase representation, and to the 0^th and 1^st order, as well as the intermediate geometrical focus of the relief-based UDS DLL model. The 0^th and 1^st order focal planes were determined by optimization for the smallest RMS Wavefront error based on the Binary 2 model, and the same positions were picked-up for analyses in the new UDS DLL model as well.

The FFT PSF results at the focal planes clearly demonstrate the difference between the order- and zone-decomposition models. While the built-in Binary 2 model only takes the selected order into consideration, the UDS DLL representation also accounts for the other out-of-focus orders creating an extensive background on the PSF. This behavior is visualized below in logarithmic scale images at the first order focus point. The False Color FFT PSF plot on the left-hand side corresponds to the order-decomposition model, while the one on the right shows the zone-decomposition results. The Cross Sections of the same PSF results at the center row, i.e. Y=0 position, are shown below.

Besides, the FFT MTF results show that the diffraction efficiency is accurately accounted for in the zone-decomposition model, but not in case of the order-decomposition.

Polychromatic results

Finally, in order to demonstrate that tracing the real ray paths in the zone-decomposition model also accounts for wavelength dispersion in the lens, we analyzed the system performance over the visible spectral range. To keep the design wavelength as a reference while extending the scope, we used the F’, e, and C’ lines by selecting the F’, e, C’ (Visible) preset. The FFT Through Focus MTF plot depicts how the first order focal plane (left hump) shifts away from the zeroth order one (hump on the right) with increasing wavelength. At the same time, at longer wavelengths more energy is diffracted to the zeroth order and less towards the first order, as expected from theory.

Summary

In this article, we demonstrated how to use User-Defined Surface DLLs to extend the capabilities of OpticStudio for modeling relief-type diffractive intraocular lenses. As IOLs are implanted close to the pupil, which serves as the STOP of the human eye, ray tracing and a one-step diffraction analysis from the exit pupil to the image plane can accurately reconstruct the light distribution diffracted into multiple orders simultaneously. We tested and verified our model by simulating an ideal bifocal IOL design in an ISO standard eye scheme. Finally, we presented the advantages of the new zone-decomposition model over the built-in order-decomposition.

References

1. A . Nemes-Czopf, D. Bercsényi, G. Erdei. Simulation of relief type diffractive lenses in ZEMAX using parametric modelling and scalar diffraction. Applied Optics, 58(32):8931-8942 (2019).
2. Ophthalmic implants—Intraocular lenses—Part 2: Optical properties and test methods, ISO 11979-2:1999.
3. A. S. Gutman, I. V. Shchesyuk, V. P. Korolkov. Optical testing of bifocal diffractive-refractive intraocular lenses using Shack-Hartmann wavefront sensor. Proceedings of SPIE - The International Society for Optical Engineering, 7718 (2010).
4. T. Eppig, K. Scholz, A. Langenbucher. Assessing the optical performance of multifocal (diffractive) intraocular lenses. Ophthalmic and Physiological Optics, 28:467–474 (2008).
5. H. Sauer, P. Chavel, G. Erdei. Diffractive optical elements in hybrid lenses: modeling and design by zone decomposition. Applied Optics, 38:6482–6486 (1999).
6. D. A. Atchison, G. Smith. Optics of the Human Eye. Butterworth-Heinemann, UK (2000).