File Exchange

image thumbnail


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


Updated 09 Mar 2021

From GitHub

View Version History

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 (2021). Optometrika (, GitHub. Retrieved .

Comments and Ratings (93)

Michele Marro

Hello Yury Petrov,

Just a question: eye model is validated?

Alexander S

Hello Lucas, Noah and Uwe,

The original code did not implement consecutive rotation. On this branch you can find a patch which implements continuos rodrigues rotations.

Regards, Alex

Alaster McDonach

Liked it and it did what I needed. Looking to use it more.

Lucas Faber

Hello Yuri,
great programm, it really helps me a lot.
I ran into a problem while applying serveral rotations around different axis to a plane and noticed that severall people had the same problem(Noah Birge,Uwe Apel).
The rotation is not visible in the .draw() function.
Is the program designed to do this? And do you have a workaround for this issue?
Myrtle42 proposed a solution on 22 Aug 2016, can I implement this without disturbing any other calculations?
Regards Lucas

Shuiqin Zheng

very very very beautiful code!! Not only useful in raytracing but also coding!

Alexander S

Good job Yury. Thanks for starting this easy to understand ray tracer.

Joseph Smalley

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?

Tanner Fadero

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.

Wilfred Mason

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

Wilfred Mason

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.

Johannes Stoerkle

Tristan Davis


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

Alaster McDonach

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?


Bradley Frank

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.

Aarón Cruz Ramírez


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!

Peter Cheimets

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

Kashmira Nakhoda

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

Stefan Schuberth

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?

Peter Cheimets

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.



Zhaozhong Chen

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?

Pratul Nijhawan

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?

Stephan Suckow

This was very helpful for me.


Olivier Ripoll

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)

Hugo Ferreira

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!

Edgar Brucke

Edvin Attebo


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!

Pascal Berto

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

Justin Little

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

Evgeniy Andreev

Evgeniy Andreev

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.


Tobias Simmet


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?


Philip Anderson

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.


Daniel Aharoni

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.

Sotirios Nousias

Sotirios Nousias


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?

Evgeniy Andreev


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!

Marshall Bremer

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.

Gerald Rothenhoefer

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').


Richard Younger


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.

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

Inspired by: hist2

Inspired: asphere parameters

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!