An independent Python framework for modeling refractive optical systems from exact vector ray tracing through wavefront analysis, optimization, and manufacturing tolerancing — built from first principles and validated against commercial tools.
Commercial optical design packages like Zemax are powerful, but they are closed boxes: you set up a system, press optimize, and trust the numbers. In addition they are limited access. I wanted to understand what happens between the prescription and the spot diagram — so I built my own optical design engine from scratch in Python, implementing every step from a ray striking a surface to a fully analyzed, optimized, and toleranced lens system.
The engine traces real, thick lenses using exact vector Snell refraction — no paraxial approximation — across field and wavelength, on general surfaces (spherical, conic, aspheric, cylindrical, flat). On top of that ray-tracing core it builds the analyses an optical designer actually uses: spot diagrams, ray and wavefront (OPD) fans, Zernike decomposition, field curvature, distortion, physical-optics PSF and MTF. It then closes the design loop with a damped-least-squares optimizer and a manufacturing-tolerancing suite, so a system can be designed, optimized, and checked for real-world yield entirely within the tool.
A guiding theme runs through the project: the connection between first-order (paraxial) optics and exact geometrical ray tracing. First-order theory tells you where an image should form and how large it should be; exact ray tracing tells you how far short of that ideal a real design falls, and why. Building the engine meant implementing both and watching them agree in the limit — which is also the clearest way to learn where the approximations break.
At its core the engine represents an optical system as an ordered list of surfaces, each carrying a radius of curvature, a thickness (the axial gap to the next surface), a glass, a clear aperture, and a shape — spherical, conic, even-aspheric, cylindrical, or flat. This prescription is the single source of truth: everything downstream, from the paraxial focal length to the final MTF, is computed by tracing light through it.
Tracing is sequential and fully three-dimensional. A ray is a position and a unit direction vector in space; a bundle of rays is launched from the object (or across the entrance pupil) and propagated surface by surface. At each surface the engine finds the ray–surface intersection by solving for where the ray meets the surface sag — a Newton iteration for curved and aspheric shapes, closed-form for planes — then applies exact vector Snell refraction using the surface normal at the intersection point. There is no paraxial small-angle approximation anywhere in the path: the full vector form of Snell’s law is evaluated, total internal reflection is detected and flagged, and rays that miss the clear aperture are marked as vignetted rather than silently dropped. The result is an honest geometric trace, valid at any field angle and aperture.
Because real glass disperses, every trace is polychromatic: the refractive index at each surface is evaluated per wavelength from the glass’s dispersion data, so chromatic effects fall out of the same machinery that produces the monochromatic spot. Alongside the geometric path, the engine accumulates optical path length — physical distance weighted by refractive index — along every ray, including the final image-space medium. That accumulated OPL, referenced to a common focus, is the wavefront, and it is what feeds the OPD maps, the Zernike decomposition, and the diffraction-based PSF and MTF.
A separate coordinate-transformation layer lets any element be tilted or decentered and then returned to the axis, which is what makes alignment tolerancing possible. And running in parallel to the exact trace is a first-order (paraxial) engine — a y–nu marginal-ray trace and the equivalent ABCD ray-transfer matrices — that provides the focal length, pupil locations, f-number, and magnification. Keeping both engines side by side is deliberate: the paraxial layer predicts the ideal, the exact trace measures the real, and the gap between them is precisely the aberration content the rest of the tool analyzes.
Once a system is traced, the engine turns the ray data into the standard analyses an optical designer uses to judge performance. Every analysis works across field and wavelength, and most can be evaluated either at each field’s own best focus or at a single fixed image plane — so the tool can report both the lens’s intrinsic quality and what a real flat sensor actually sees.
Geometric image quality. Spot diagrams show the footprint of a traced ray bundle at the image, overlaid with the diffraction (Airy) disk for scale, and condensed into an RMS spot radius versus field. Ray fans plot transverse ray error across the pupil, making the signature of each aberration directly visible in the shape of the curve.
Wavefront analysis. From the accumulated optical path, the engine produces exit-pupil OPD maps and the RMS wavefront error in waves, referenced to the Maréchal (diffraction-limited) criterion. A Zernike decomposition breaks the wavefront into named aberrations — defocus, astigmatism, coma, spherical — so their growth with field angle can be tracked individually.
Field-dependent behavior. Field curvature is measured as the shift of best focus with field angle; distortion as the percentage departure of chief-ray image height from the ideal; telecentricity as the chief-ray angle in object and image space. Together these characterize how a design behaves away from the axis and on a flat detector.
Diffraction and image quality. A physical-optics layer computes the point-spread function and Strehl ratio, the modulation transfer function (MTF), and encircled energy — the metrics that actually decide resolution. Because these accept an explicit evaluation plane, the MTF can be reported at the flat sensor including field-curvature defocus, not just at each field’s private best focus.
Chromatic analysis. Axial (longitudinal) and lateral (transverse) color quantify how focus and image height drift with wavelength, and chromatic focal shift traces focus across the spectrum — the measures that drive glass selection and achromatization.
Energy, coatings, and polarization. Because each ray carries a flux-normalized electromagnetic field rather than only a geometric path, the engine also models the physical layers that set real throughput and contrast. Thin-film coatings are evaluated with the Fresnel equations and a transfer-matrix method, giving per-surface reflectance and transmittance versus wavelength and angle — including anti-reflection stacks. Polarization ray tracing propagates the Jones vector through the system, reporting diattenuation and retardance across the pupil and a polarization-resolved (vector) PSF.
And ghost-reflection analysis enumerates double-bounce specular reflection paths — the brightest contributors to veiling flare — and ranks them by irradiance, so the benefit of an anti-reflection coating can be quantified directly as a suppression factor. (This is the ghost component of stray light, not full non-sequential stray-light modeling.)
A design is only useful if it survives fabrication and assembly, so the engine extends the optimization loop into manufacturing. Every construction parameter can be perturbed within a realistic tolerance — surface radius, element thickness, and refractive index on the fabrication side; element tilt, decenter, and despace on the alignment side. A sensitivity analysis perturbs each one in turn and ranks them by their effect on performance (RMS wavefront, spot, boresight, or divergence), identifying which few tolerances actually drive the design. A Monte-Carlo analysis then randomizes all of them together over hundreds of virtual builds and reports the yield — the fraction of assembled systems that still meet a chosen specification. A compensator — an adjustment a technician would tune at assembly — can be folded in: the built-in back-focus refocus, or, for the afocal expander, a re-collimation gap adjustment, so the predicted yield reflects the real build process. In the beam-expander study, modeling a single gap compensator is what lifts the as-built yield to an acceptable level, exactly as it would on the bench.
A ray tracer that only analyzes a fixed prescription is an evaluation tool; to be a design tool it has to close the loop and improve the prescription automatically. The engine does this by exposing every constructional parameter — surface curvatures, element thicknesses, air gaps, stop position — as an optimization variable, and letting a numerical optimizer drive them toward a target.
The design workflow mirrors how a lens is actually developed. A starting system is seeded from first-order theory — focal lengths and power distribution from the paraxial equations, radii from thin-lens formulas, spacings from the layout — then traced and analyzed to expose its aberrations. The optimizer then refines it: the designer declares which parameters are free, which performance quantities are targets (held at a value) and which are objectives (minimized), and the optimizer searches the variable space for the best compromise. The process iterates — analyze, adjust the setup, re-optimize — exactly as in commercial practice.
Two optimizers are provided. A scalar-metric optimizer minimizes a single weighted figure of merit, and a damped-least-squares engine minimizes a vector of operands directly — the method commercial optical design software uses. The merit is multi-operand and mixes constraints with objectives: first-order quantities such as focal length, f-number, and magnification can be held while geometric spot size, RMS wavefront error, field curvature, distortion, telecentricity, or beam divergence are minimized, each with its own weight, target, or inequality bound, and each variable can be bounded to a physical range. Critically, the ray-based operands are plane-aware — they can be evaluated at a fixed flat sensor rather than at each field’s private best focus — so the optimizer corrects for what a real detector sees, including field-curvature defocus.
In practice this drives real design decisions. For the camera triplet, holding EFL and f/number while minimizing wavefront and field curvature — and adding a field-flattener element between optimization passes — pulled the flat-sensor corner spot from 80 µm to 16 µm. For the beam expander, minimizing output divergence while holding the 5× magnification collimated the beam from 2.6 mrad to 0.06 mrad. In each case the optimizer did not invent the design; it refined a physically-seeded starting point toward a specified, measurable goal.
A custom solver is only trustworthy if its outputs can be checked against something independent. The engine was validated at three levels — internal consistency, analytical benchmarks, and a direct comparison with commercial design software.
First-order consistency. The exact 3D ray trace must reproduce paraxial predictions in the limit where paraxial theory is exact — small apertures and fields. Effective focal length, f-number, pupil positions, and magnification computed from the real ray trace are checked against an independent y–nu paraxial engine and the equivalent ABCD matrices, and agree to numerical precision. This confirms the ray–surface intersection and refraction core before any aberration is interpreted.
Analytical benchmarks. Individual results are compared to closed-form optics: the diffraction spot matches the Airy radius 1.22·λ·(f/#); the RMS-wavefront-to-Strehl relationship follows the Maréchal criterion; the achromatic doublet, designed to the Σφ/V = 0 condition, shows the predicted collapse of axial color relative to an equivalent singlet; and the beam expander’s residual divergence bottoms out at the Gaussian far-field floor λ/(πw), never below it.
Comparison with commercial software (Zemax OpticStudio). The strongest check reproduces a published, independently computed design: Edmund Optics’ stock-optic 5× Galilean beam expander (−12 mm and 60 mm catalog lenses, 632.8 nm HeNe, 3 mm input beam), optimized — as in their Zemax example — by adjusting only the lens separation. Tracing the identical prescription, the engine recovers the expansion, the collimated output, and the exit-pupil wavefront to within ~3% of Zemax’s reported values:
Metrix, Zemax Optic Studio, This Engine, Difference: Wavefront P-V (Waves), 0.4630, 0.4740, 2.4% ; Wavefront RMS (waves), 0.1311, 0.1360, 3.7%
Agreement at this level — a from-scratch tracer matching OpticStudio’s peak-to-valley and RMS wavefront on the same design — demonstrates that the ray tracing, optical-path accumulation, and afocal wavefront referencing are all implemented correctly end to end.
The engine was developed alongside a series of real design problems, each chosen to stress a different part of the tool. Rather than one showcase lens, the case studies demonstrate that the same framework handles imaging and afocal systems, monochromatic and polychromatic design, and on-axis and strongly field-dependent performance.
Achromatic doublet. A cemented crown–flint doublet designed to the Σφ/V = 0 achromatic condition, compared directly against an equivalent single lens. It exercised the polychromatic trace and the axial-color analysis, and confirmed the predicted collapse of chromatic focal shift relative to the singlet — the engine’s chromatic machinery validated against textbook theory.
Full-frame camera triplet. A fixed-focal-length f/8 lens for a 36 × 24 mm sensor over a 40° field. This was the field-dependent stress test: it drove the Zernike-vs-field analysis, the flat-sensor evaluation, the two-stage optimization workflow, and the addition of a Petzval field-flattener — reducing the flat-sensor corner spot from 80 µm to 16 µm and exposing field curvature as the true limit of a bare triplet.
5× beam expander (1550 nm). A Galilean afocal expander that widens a 1550 nm Gaussian beam fivefold within a 75 mm track. It exercised the afocal image-space mode (angular spot in mrad, plane-referenced wavefront) and the divergence operand: optimizing the lens radii and separation while holding the 5× ratio collimated the collimation output from 2.6 to 0.06 mrad. It also drove the full tolerancing stack — a sensitivity ranking of alignment and fabrication tolerances, and a Monte-Carlo yield that rose from 34% to 99% once a single gap-spacing compensator was modeled.
import optray as ot
lens = ot.System([486.1, 587.6, 656.3], object_distance=float('inf'), field_type='angle')
lens.set_primary(1)
lens.add(radius=72.05, thickness=3.60, glass='N-BK7', semi_dia=3.75, stop=True)
lens.add(radius=-21.90, thickness=1.10, glass='F2', semi_dia=3.75) # cemented flint
lens.add(radius=-57.01, thickness=75.0, glass='air', semi_dia=3.75)
lens.set_fields([(0, 0)]); lens.focus()
res = lens.trace(ot.Pupil(n_rays=2000))[0]
res.axial_color() # per-wavelength focus shift (span 38 µm vs 1152 µm singlet)
res.polychromatic_spot() # white-light spot + Airy diskCamera triplet — damped-least-squares optimization at the flat sensor
lens.optimize_dls(
vary=['s1.radius', 's2.radius', 's5.radius', 's6.radius', 's6.thickness'],
targets={'efl': 49.45, 'fnumber': 8.0, # held
'wavefront': 'min', 'field_curvature': 'min'}, # minimized at plane='sensor'
plane='sensor', source=ot.Pupil(n_rays=400))
res = lens.trace(ot.Pupil(n_rays=2000))
res[2].primary.mtf(plane=res[0].primary.best_focus) # corner MTF on the flat sensorTolerancing — sensitivity ranking + Monte-Carlo yield with a compensator
ELEMS = [(1, 2), (3, 4)]
FAB = {'radius': 0.1, 'thickness': 0.05, 'index': 0.0005}
lens.sensitivity(ELEMS, fab_tolerances=FAB).tornado('rms_waves') # which tolerance hurts most
mc = lens.montecarlo(ELEMS, fab_tolerances=FAB, n_trials=300, seed=0, compensate='focus')
print(mc.yield_fraction('rms_waves', 0.07)) # as-built yield against the specThe engine is written in pure Python on a scientific stack of NumPy, SciPy, and Matplotlib, with no dependency on any optical-design software — every result comes from first-principles code. It is about 5,400 lines across 17 modules (ray-trace core, surfaces, analyses, physical optics, optimization, tolerancing, coatings, polarization, ghosts), and is backed by a regression suite of roughly 230 tests that check the physics against analytical and paraxial references.
Its scope is deliberately bounded. The engine is a sequential ray tracer: light propagates through surfaces in a defined order, which covers the vast majority of lens and imaging-system design — refractive and reflective elements, spherical, conic, aspheric, cylindrical, and flat surfaces, across field and wavelength. It is not a non-sequential tool, so general stray-light and scatter analysis (housing reflections, surface-roughness scatter, baffles) is outside its scope; the one exception is double-bounce ghost-reflection analysis, which is handled explicitly. It also does not model diffractive or gradient-index surfaces, or thermal and mechanical effects.
Python / NumPy · ~5,400 lines, 17 modules, ~230 physics tests · sequential ray trace → wavefront, PSF/MTF, optimization & tolerancing · validated against Zemax OpticStudio to ~3%.
Selected implementations are presented above. The complete source code, including the ray-tracing core, optical analysis, optimization algorithms, and validation tests, is maintained in a private GitHub repository