File Exchange

image thumbnail


version (375 KB) by Yury
Optometrika MATLAB library implements analytical and iterative ray tracing approximation to optical image formation using Snell’s and Fresne


Updated 23 Jun 2020

GitHub view license on GitHub

Developer: Yury Petrov

Optometrika library implements analytical and iterative ray tracing approximation to optical image formation using Snell’s and Fresnel’s laws of refraction and reflection.
Currently, the library implements refractive and reflective general surfaces, aspheric (conic) surfaces with astigmatism, Fresnel surfaces, cones and cylinders (elliptic too), planes, circular and ring-shaped apertures, rectangular flat screens, spheroidal screens, and a realistic model of the human eye with accommodating lens and spheroidal retina. See example*.m files for examples of ray tracing in general (user-defined shape) lenses, aspheric lenses, Fresnel lenses, prisms, mirrors, and human eye.

The library traces refracted rays, including intensity loss at the refractive surface. Reflected rays are currently traced for mirrors and also for a single total internal reflexion or double refraction, if it happens. Note that the Bench class object is not a real physical bench, it is only an ordered array of optical elements, and it is your responsibility to arrange optical objects in the right order. In particular, if you need to trace rays passing through the same object multiple times, you have to add the object multiple times to the bench array in the order the object is encountered by the rays. For example, double refraction/reflection for cylindrical and conical surfaces can be calculated by adding the surface twice to the bench.

The library is very compact and fast. It was written using Matlab classes and is fully vectorized. It takes about 2 seconds to trace 100,000 rays through an external lens and the human eye (8 optical surfaces) on a 3 GHz Intel Core i7 desktop. Fresnel lens tracing is somewhat slower due to looping through the Fresnel cones describing the lens surface. Tracing through user-defined (general) surfaces is significantly slower due to iterative search of ray intersections with the surface.

Thank you for downloading Optometrika, enjoy it!

List of examples:

example1.m: tests the basic functionality of the Optometrika library
example2.m: demonstrates the Optometrika's optical model of the human eye
example3.m: demonstrates accommodation of the human eye by minimizing the retinal image
example4.m: tests a ring lens with the cosine surface profile defined in coslens.m
example5.m: tests planar mirrors
example6.m: tests planar and parabolic mirrors (a Newtonian refractor telescope)
example7.m: tests a Fresnel lens
example8.m: tests a lens with polynomial aspheric terms
example9.m: tests cone mirrors
example10.m: tests cylinder and cone surfaces with double refraction
example11.m: demonstrates ray tracing for rays originating inside the human eye
example12.m: draws a lens and determines its front surface, back surface, and total height. Makes an animated gif of the lens and an engineering drawing of the lens.
example13.m: tests refraction through the lens edge and backward rays refraction (sub-aperture Maksutov-Cassegrain telescope)
example14.m: tests refraction through a lens with astigmatism (different vertical and horizontal radii of curvature)
example15.m: simulates a hexagonal array of spherical micro lenses
example16.m: demonstrates STL export of various lenses

Cite As

Yury (2020). Optometrika (, GitHub. Retrieved .

Comments and Ratings (87)

Hi Yury,

When you say it takes 2 seconds to trace 100,000 rays through 8 surfaces, are you referring to 100,000 source rays, or a total of 100,000 ray segments? Approximately how many times slower is ray tracing through user-defined surfaces?

Hi Yury,

I love using Optometrika, thanks so much for putting so much time into writing it!

I work with light sheet fluorescence microscopy, where I often need to ray trace through plano-convex cylindrical lenses. Your CylinderLens.m classdef is very useful, but I was wondering if there was a way to introduce an argument that would "crop" (i.e. introduce a bounding box or maximum outer diameter) the CylinderLens.m object in the bench container, similar to how you can specify an "outer diameter" argument for the Lens.m object. This would help me with figure generation, since the cylindrical lens wouldn't be "intersecting" the rays at different points in the optical path (even though they don't refract through the lens at that point).

Thanks again!


Wilfred Mason, Optometrika is not my main line of work, so I didn't get to writing a tutorial. But I feel that if you understand how Matlab works in general (e.g., understand its data structures an know how to debug a program) you should find it fairly easy to understand Optometrika by following the provided examples. Any time you are in a doubt about what is going on just put a breakpoint and dive into the Optometrika function that you want to understand. Otherwise use "help function_name" to read its description. I tried to put enough information there.

Can be very useful but the lack of documentation makes it difficult to use. Please consider incorporating a documentation package.

Awesome package! Does anyone know how to modify the input rays so that they have a Gaussian distribution. I see the functions, but I don't know how to incorporate them.


I am having trouble getting this to work with freemat? Any chance someone knows which parts of these would hold up the code? It has to be formatting of some kind. In a general sense, is there specific items that people have ran into with this?

Nick Stuart

Enjoyed using it.

Jin Yang

Nice job! I was wondering is that easy to modify your codes to trace light rays pass through non-spherical (parametric function is given) transparent lens ?

Noah Birge

Hi Yury,

I'm having some trouble making my aperture rotate.

a=Aperture(positionVector, rectangularDimensions);
a.rotate([0 1 0],pi/2);

^ when I use "a.display()" in the command line, I get the correct information, but bench.draw() shows the aperture in the un-rotated position. The rays seem to act as though the aperture is in the correct location, but the visualization of the aperture location seems to be different. Am I doing something incorrectly?


Stefi Angu

Dear Yury, I have a question on how to get number of rays in a desired coordinates/screen to represent the light intensity.Thank you.


Stefi, reflections are implemented, see examples 5, 6, 9, 13. You can build a retroreflector, e.g., a corner reflector or a refractive sphere using the existing objects like Planes arranged into a 3D corner.

Stefi Angu

Dear Yury, I'm trying to simulate a confocal system which requires Retro-reflection function in Optometrika. Is there anyway i can capture the reversed reflection because i cant seem to make it work.


f d, export_stl() works with Plane surfaces in Optometrika 2.2. Also, the internal diameter (hole) was added to Plane. Thanks for pointing this out!


Konstantina, beamsplitters are not implemented. A possible workaround is to run the same beam of rays through a bench twice: the first time with a mirror, the second time - without. Then you can combine the two sets of Rays using Rays.append() function and send the resulting bundle to a screen or elsewhere. Just note that no phase information is currently taken into account in Optometrika, it is a purely geometric optics approximation.


Uwe Apel, you can use the Retina object as a spherical screen.


Wick, thanks a lot for catching this bad bug! This problem was fixed in Optometrika 2.2.


Peter Cheimets, the bug affecting backward rays and hyperbolic surfaces was fixed in Optometrika 2.2. Thanks for finding it!

I'm trying to simulate a Cassegrain telescope. I started with example 13 (Maksutov-Cassegrain), removed the lenses and tweaked the element positions so they focus. Then I put in a -1 conic constant for both mirrors. This worked fine. However, in a Cassegrain, the secondary has a <-1 conic constant. When I changed the conic constant to less than -1 (even as small a change as going to -1.00000001), the beams coming off the secondary no longer converged. Am I missing something, or is this incapable of dealing with conic constants less than -1? (I downloaded my version in April.

Thanks Peter

Hi Yuri,
Thank you for developing this toolbox, it is extremely powerful.
I am trying to use the Plane element in your package to transform it as a coded aperture of 0 and 1 and add a source that is an extended object.
Have you ever tried introducing such elements before?


Hi Yury! Are beamsplitters implemented in the optical elements and if not, can they be implemented? Thank you for your answer! Konstantina

Is there a solution to Wick's problem?


I have found a bug in using an aspheric lens with rays.intersection. Consider this snippet of code:


R = 50; k = -1.3; a4 = 5e-7;

bench = Bench;
bench.append(AsphericLens([0 0 0],50,R,k,[0 a4],{ 'air' 'n-bk7'}));
bench.append(Screen( [20 0 0 ], 50, 50, 512, 512 ));

rays_in = Rays( 2, 'collimated', [ -5 0 0 ], [ 1 0 0 ], 2*8.789, 'linear', 'air', 532e-9); % just two rays in the XY plane at +/- 8.789 above the centerline
rays_out = bench.trace( rays_in );
pass_angle_relative_to_centerline = atand((rays_out(3).r(:,2) - rays_out(2).r(:,2)) ./ (rays_out(3).r(:,1) - rays_out(2).r(:,1))) % arctan (deltaY / deltaX)

The angle relative to the centerline after passing into the lens surface at +/- 8.789 is -/+ 6.4 degrees. However, at that height, a horizontal ray has an angle of incidence of 10 degrees. For n = 1.519 glass, the transmission angle relative to the internal surface normal is 6.56 degrees. So the transmitted angle relative to the centerline should be -/+ 3.44 degrees.

The points where the rays intersect the aspheric surface appear to be correct so the problem must be with the calculation of the angle.

I see you had a check for the Optimization toolbox in that section of the code. Disabling the Optimization toolbox did not change the result. Thank you in advance for taking another look at your aspheric surface calculations.


Thanks for the submission. I'm just now getting started. Is there any built-in feature to determine wavefront error (optical path difference, OPD) shy of me evaluating the optical path length for each ray? Thanks.

f d

Dear Yury, I encounter a problem, the front lens is aspheric, and the back len is planar, when I use export_stl() to export, cause error below,
error: export_stl (line 79)
p2 = asphlens( Y2 * cs, Y2 * sn, pars2, 0 );
no variable pars2.

I found is the problem of planar back len whose D is scalar not vector, how to solve this question?

I downloaded the program last week and have been playing with it. It looks very powerful. Has anyone written up some notes on how to use it? I have been digging into the software itself, and looking at the examples, but I was hoping for something more manual-like.



Uwe Apel

Hi Yury - After switching from Matlab2013b to Matlab2018 now almost all works fine. However, I would like to use a spherical screen, but I wasn't able to find any example how the default setting R=inf could be changed successfully. Could you kindly advise? Many thanks.

Ann Oman

Hi Yury - What are some of the anticipated use cases for the eye model? We are curious to test focusing on near and far objects and what is seen by the eye as well as potential near and farsighted cases.

Is this possible?

Thanks a lot for such a wonderful library and so much respect for the amount of help you did to the society. Is the 3D tracing of rays through an array of some points, say 3*3 matrix of points in a 3D space and obtaining the final outputs as zones on image screen possible with this library, say a set of some rays are passing point by point to lens's 3D surface, then the obtained image on screen would be in form of zones, so I needed some suggestions regarding that. In the designed library, which examples I should refer for the same?

Kind regards!


Dear 123, run "help coslens" to see a description of how to define a freeform surface to be used with Optometrika. It doesn't have to be described by a formula, like in example4, you just need a function that provides the lens surface's height and normal at any given (x,y) point in the lens plane. It could be obtained by spline interpolation of the surface data from a look-up table, for example.


Yu Xiao

Zhijie Ma


Hello Yury,
This program is very usefull. Can this program use for freeform lens with the coordinates (X, Y, Z) that is produced by the user? Not just use for the 'GeneralLens' 'AsphericLens' 'FresnelLens', etc. I have sent an email to you in this web from the email Contact. I hope you will see it.
Thank you very much.



Hello Yury,

There is bug in example16.m when use matlab 2018b. It does not work.

When the example is used in freeform lens which only has coordinate information, how to modify it? In other words, how to use them to simulate freeform lens and reflectors?

Thank you very much.

Uwe Apel

Hello Yury,
Just another question: All example solutions that use function "quadric_intersection" do not work, failure is mismatch of matrix dimensions at both last formulas:
t1 = 0.5 * ( -b - sqrt( D ) ) ./ a;
t2 = 0.5 * ( -b + sqrt( D ) ) ./ a;

r1 = r0 + n .* t1;
r2 = r0 + n .* t2;
t1, t2 are only one-dim vectors while r0 and n are 127x3 matrices. I'm working with Matlab13b,
might this be a version issue? Any hint? Thanks.

Uwe Apel

Helly Yury,
Very useful approach for visualization of optical test benches. Some more basic information would be appreciated.
I observed some not-understood behaviour. So I modified example5 by rotating the 2nd mirror additionally a small amount around the y-axis. The course of the rays behaves in plausible way, but the orientation of the displayed mirror plane is completely different:
mirror2 = Plane( [ 60 50 0 ], 40, 40, { 'mirror' 'air' } ); % pay attention to the glass order here!
mirror2.rotate( [ 0 0 1 ], -pi / 4 );
mirror2.rotate( [ 0 1 0 ], -pi / 18 );
bench.append( mirror2 );

Any idea?

Jongmin Kim

This is very good example!!

Joshua Bone

What is the purpose of the attenuation property of the Ray object?

This was very helpful for me.


you may want to replace the for loop code adding the polynomials in the aspheric.m by a call to polyval, or rewrite it using Horner's algorithm to make it faster and more reliable (polyval uses the Horner algorithm)

Hello Yury. First of all, great work with this toolbox! I read the comments here and it seems there was a grating simulation tool in a previous version. Why was it removed, and can you provide me with that tool? Thanks!


Hey Yury - thanks for the quick reply. I'm on R2016a so that might not be the issue, but I'll keep an eye out for your next update. Thanks!


If you have problems similar to the one described by Justin, R S, and Pascal, please update your Matlab to 2016 or later. I'll fix this in the next update, so it will work for earlier Matlab versions.


I have the same problem as Pascal & Justin, any help would be appreciated. Thanks!

Hello Yury,
I have the same problem than Justin Little...
Please, do you know where it could come from?

Hello Yury,

I get the following error when attempting to trace rays in example1 (and other examples). I am using Matlab 2015B. Any thoughts?

Error using ./
Matrix dimensions must agree.

Error in Rays/intersection (line 712)
nrms = nrms ./ sqrt( sum( nrms.^2, 2 ) );

Error in Rays/interaction (line 732)
[ rays_out, nrms ] = self.intersection( surf );

Error in Bench/trace (line 327)
rays( i + 1 ) = rays( i ).interaction( self.elem{ i }, out_fl );

Error in example1 (line 42)
rays_through = bench.trace( rays_in ); % trace rays

Hello, Yury

First I want to say thank you for the great toolbox! It really helped me in studying and scientific work.
What do you think about adding the project to GitHub? I'd like to contribute to it, like many others I believe.



Hello Yury,

I am trying to add an additional rotation to mirror1 in example5:
(ZIP file downloaded on 9/21/2017)

% planar mirror
mirror1 = Plane( [ 60 0 0 ], 40, 40, { 'air' 'mirror' } ); % pay attention to the glass order here, it defines the mirror orientation!
mirror1.rotate( [ 0 0 1 ], -pi / 4 );
mirror1.rotate( [ 0 1 0 ], 0.00000001 ); % add this line to example5.m
bench.append( mirror1 );

The rendered image draws mirror1 at an unexpected angle.

What is the correct way to add additional surface rotations?



Marie Lu

Hi Yury,
great software! What is the default unit for lengths ?
Thanks a lot!

Hello everyone,
I am having a bit of trouble trying to grasp how to use the toolbox. Is there a simple user manual that I can use to better understand the toolbox?


Hello Yury,

The code Just what I was looking for to get solid angle sensitivity of a photosensor buried in clear plastic.

One question: I wish to trace rays that go through AND past a hemisphere lens, so I have been using
bench.trace( rays_in, 0 );
to allow missed rays t be traced. Works fine if I use a "conelens" (but this is not what I need) , but not with "lens" (which I use to build the 1/2 sphere.
Looking in Rays, the line
rays_out.I( in, : ) = self.I( in );
or its equivalent is not included for "lens", (as it is for "conelens" and others) so the intensity is not passed on.
Is there a fix for this please (not trivial because there is no obvious indices 'in', Otherwise the code is very close to ideal for me.


Excellent optics package! I think there is an error in the focal point calculation in Rays.m. I think f = osr + +d * nav; should actually be f = osr + rav +d * nav;

Noam Cohen

Could be even greater if will be expanded to NIR/SWIR/MWIR/LWIR.


Great toolbox! Thank you, Yury!

Question for anyone; does anyone know how to add a custom lens like a Powell lens? I was able to generate one with a different script, but the rays don't interact with it when I incorporate it into an optical system. Thanks in advance!


In case anyone runs into the same problem, I corrected it by editing the rotate method in Surface.m:

Comment out lines:

Add in:
% Calculate new normal vector to the plane
self.n = rodrigues_rot(self.n, rot_axis, rot_angle);

% Calculate rotation between original and final vectors
rotvec = vrrotvec([1 0 0],self.n);

% Assign correct rotation axis and angle
self.rotax = rotvec(1:3);
self.rotang = rotvec(4);

I didn't check lines 44-49 are still correct since I don't rotate more than 180 deg

Clarification welcome if there is a better way to apply multiple rotations in 3d


How to create a prism? Can't find it in the examples. Thanks!

Max Lowndes

I'm not really clear on how to use this! Would love some help as it seems to do exactly what I need. Trying to run a ray tracing simulation of a fiber optic!

xizhi gu

Hi,Yury, Thanks for the toolbox!
I have a question that why the "grating" is removed from the newest version? Because I can't find it anymore..

Robert Shah

Whats the default length unit, mm?


Yury, your work is impressive. Is it possible to treat influence on the image of diffraction grating placed between object and imaging lens?


Very impressive and useful, thanks for sharing!

Thank you. This is exceptional.


Yury, Thanks for sharing this AWESOME tool.
I was sitting down to write something for doing the same job... and would not have done it nearly this well. A nice education in OOP.

Juan Perez

Hi, this is a very useful and versatile tool, congratulations and thanks for sharing it! I am trying to get the coordinates of the intersection of the rays with the retina (in example2.m for example), can you give me some idea of how to do this?


Awesome ideas and code. Documentation could be much less sparse.


If you want to use a material which is not in the refrind.m file the easiest is to edit the file. Just add another 'case', i.e., case 'quartz' nref = [ 1.4585 1.4631 1.4564 ];
Just make sure that the order of the refractive indices corresponds to wavelengths of 587.6, 486.1, and 656.3 nm.

Hi, my first impression is very good. Great work!
But, it takes some round about to figure out which glass materials are available and which properties are assumed.
The links given under doc(refrindx) seem to be dead.
So, the only way to get the information is by edit(refrindx).
Would be nice to make the refractive index an input argument to the lens function (as an alternative to 'glass').



Excellent program. I am wondering if there is a solution to create two source rays in one model. Tried to fix this by using two aperatures, however they block each other. Keep up the good work.


Excellent. It would be nice to implement cylindrical lenses


Excellent and makes wonderful use of OOP, using a virtual optical bench object.


Fixed a bug in Aperture.m incorrectly displaying rectangular aperture orientation on the bench

Fixed a bug preventing displaying a point image on a screen

Bug fixing

Fixed incorrect surface normals calculation for surfaces with polynomial terms, added inner diameter (a hole) to the Plane surface , extended export_stl() to Plane surfaces.

Fixed a bug affecting backward propagating rays interacting with hyperbolic surfaces.

Maintained as a GitHub repository now.

Changed export_stl() function to work correctly with the stlwrite() function bulid into Matlab 2018b.

Updated the export_stl() function to export GeneralLens class surfaces.

Added export_stl() function for lens surface export in STL format. See example16.
Various bug fixes

Various bug fixes

Added several new materials
Fixed several bugs, including the one preventing the library working on pre-2016 Matlab versions

Added STL and DXF export of a lens's 3D shape and 2D profile respectively to the draw_lens_engineering() function. Install stlwrite() and DXFLib to use this functionality.

Minor bug fixes.

Added engineering drawing of Fresnel lenses to draw_lens_engineering(). See example7 for an example.
Fixed a bug in Fresnel surface rotation.

Added rectangular apertures. See example4.
Fixed various bugs in the FresnelLens implementation

Implemented astigmatic lens surfaces: use a 2x1 vector for the radius of curvature (see example14).
Implemented elliptical cone/cylinder surfaces: use a 2x1 vector for the cone/cylinder diameter (see example10).
Various bug fixes

Implemented ray tracing for backward-propagating rays (example13)
Various bug fixes

Bug fixes

Added flange option to the draw_lens_engineering() function.

Various bug fixes.

Added draw_lens_engineering() function to automatically generate engineering drawing for lens manufacturers. Currently only lenses without holes are implemented. See example12 for a demo.

Example12 illustrates the use of lens_dims() function to draw a lens (full surface) and calculate its sags and volume.
Various bug fixes.

For a new way to display rays, see example11.
Added density information for many materials in refrindx.m. Use lens_dims() function to get lens volume.
Example11 demonstrates tracing from inside the human eye.
Various bug fixes.

Total internal reflection enabled. Various bug fixes related to multiple refraction/reflexion.

Added cylindrical surfaces, including double refraction on such surfaces (example 10). Fixed some bugs.

Added a cone-shaped surface ConeLens()
Added ring-aperture option for all surfaces except the FresnelLens.
Rays missing optical elements are not discarded any more but traced to the next element.
General bug fixes

Added AsphericLens class describing strongly aspheric surfaces (conic + even polynomial terms, example8).
Added a fast focal_point() function looking for the bundle's focal point without iterations (example1).

Added a Fresnel lens implementation. Added operations on a bench, such as bench rotation, element prepending and removal.

Replaced the missing Bench.m file.

Various bug fixes and improvements. Spheroidal screens are implemented now by the Retina object.

Some bug fixes and improvements, added spheroidal screens implemented by the Retina class.

Added basic documentation to the classes.

Added: (i) appending ray bundles to an existing bundle, (ii) mirror elements (general, conic, planar). To specify a surface as a mirror use 'mirror' glass specifier (see example5.m and example6.m).

Small bug fixes.

Added tracing through general (user-defined) surfaces (see example4.m).

Fixed a bug affecting ray tracing through strongly hyperbolic surfaces.

Fixed a bug which affected ray tracing through extremely hyperbolic surfaces.

MATLAB Release Compatibility
Created with R2018b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Inspired by: hist2

Inspired: asphere parameters