Wavefront Error Optimization

import numpy as np

from optiland import optic, optimization

Define a starting lens:

lens = optic.Optic()

# add surfaces
lens.surfaces.add(index=0, radius=np.inf, thickness=np.inf)
lens.surfaces.add(index=1, radius=50, thickness=3, material="SK16")
lens.surfaces.add(index=2, radius=-100, thickness=5)
lens.surfaces.add(index=3, radius=-100, thickness=3, material=("F2", "schott"))
lens.surfaces.add(index=4, radius=50, thickness=5, is_stop=True)
lens.surfaces.add(index=5, radius=50, thickness=3, material="SK16")
lens.surfaces.add(index=6, radius=-20, thickness=40)
lens.surfaces.add(index=7)

# set aperture
lens.set_aperture(aperture_type="EPD", value=10)

# set fields
lens.fields.set_type(field_type="angle")
lens.fields.add(y=0)

# set wavelengths
lens.wavelengths.add(value=0.48)
lens.wavelengths.add(value=0.55, is_primary=True)
lens.wavelengths.add(value=0.65)

lens.draw()

Define optimization problem:

problem = optimization.OptimizationProblem()

Add operands (targets for optimization):

"""
Add a wavefront error operand for all wavelengths.

Use Gaussian quadrature distribution for the rays (see distribution documentation for
more information).
"""

for wave in lens.wavelengths.get_wavelengths():
    input_data = {
        "optic": lens,
        "Hx": 0,
        "Hy": 0,
        "num_rays": 3,
        "wavelength": wave,
        "distribution": "gaussian_quad",
    }
    problem.add_operand(
        operand_type="OPD_difference",
        target=0,
        weight=1,
        input_data=input_data,
    )

Define variables - let all radii vary (surface indices 1 through 6)

for surface_number in range(1, 7):
    problem.add_variable(lens, "radius", surface_number=surface_number)

Check initial merit function value and system properties:

problem.info()

╒════╤════════════════════════╤═══════════════════╕
│    │   Merit Function Value │   Improvement (%) │
╞════╪════════════════════════╪═══════════════════╡
│  0 │                  24631 │                 0 │
╘════╧════════════════════════╧═══════════════════╛
╒════╤════════════════╤══════════╤══════════╤══════════╤══════════╤════════════════════╕
│    │ Operand Type   │   Target │   Weight │    Value │    Delta │   Contribution (%) │
╞════╪════════════════╪══════════╪══════════╪══════════╪══════════╪════════════════════╡
│  0 │ OPD difference │        0 │        1 │ 104.534  │ 104.534  │            44.3647 │
│  1 │ OPD difference │        0 │        1 │  89.8916 │  89.8916 │            32.8063 │
│  2 │ OPD difference │        0 │        1 │  74.9868 │  74.9868 │            22.8291 │
╘════╧════════════════╧══════════╧══════════╧══════════╧══════════╧════════════════════╛
╒════╤═════════════════╤═══════════╤═════════╤══════════════╤══════════════╕
│    │ Variable Type   │   Surface │   Value │ Min. Bound   │ Max. Bound   │
╞════╪═════════════════╪═══════════╪═════════╪══════════════╪══════════════╡
│  0 │ radius          │         1 │      50 │              │              │
│  1 │ radius          │         2 │    -100 │              │              │
│  2 │ radius          │         3 │    -100 │              │              │
│  3 │ radius          │         4 │      50 │              │              │
│  4 │ radius          │         5 │      50 │              │              │
│  5 │ radius          │         6 │     -20 │              │              │
╘════╧═════════════════╧═══════════╧═════════╧══════════════╧══════════════╛

Define optimizer:

optimizer = optimization.OptimizerGeneric(problem)

Run optimization:

optimizer.optimize()

  message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
  success: True
   status: 0
      fun: 0.4929267697216927
        x: [-5.850e-01 -1.977e+00 -2.011e+00 -4.529e-01 -4.369e-01
            -2.186e+00]
      nit: 9
      jac: [ 1.326e-01  9.691e-02 -2.440e-01  4.657e-01  9.640e-02
             4.481e-02]
     nfev: 70
     njev: 10
 hess_inv: <6x6 LbfgsInvHessProduct with dtype=float64>

Print merit function value and system properties after optimization:

problem.info()

╒════╤════════════════════════╤═══════════════════╕
│    │   Merit Function Value │   Improvement (%) │
╞════╪════════════════════════╪═══════════════════╡
│  0 │               0.492927 │            99.998 │
╘════╧════════════════════════╧═══════════════════╛
╒════╤════════════════╤══════════╤══════════╤══════════╤══════════╤════════════════════╕
│    │ Operand Type   │   Target │   Weight │    Value │    Delta │   Contribution (%) │
╞════╪════════════════╪══════════╪══════════╪══════════╪══════════╪════════════════════╡
│  0 │ OPD difference │        0 │        1 │ 0.485152 │ 0.485152 │           47.7499  │
│  1 │ OPD difference │        0 │        1 │ 0.155018 │ 0.155018 │            4.87508 │
│  2 │ OPD difference │        0 │        1 │ 0.483243 │ 0.483243 │           47.375   │
╘════╧════════════════╧══════════╧══════════╧══════════╧══════════╧════════════════════╛
╒════╤═════════════════╤═══════════╤═══════════╤══════════════╤══════════════╕
│    │ Variable Type   │   Surface │     Value │ Min. Bound   │ Max. Bound   │
╞════╪═════════════════╪═══════════╪═══════════╪══════════════╪══════════════╡
│  0 │ radius          │         1 │   41.5009 │              │              │
│  1 │ radius          │         2 │  -97.652  │              │              │
│  2 │ radius          │         3 │ -101.134  │              │              │
│  3 │ radius          │         4 │   54.7111 │              │              │
│  4 │ radius          │         5 │   56.3141 │              │              │
│  5 │ radius          │         6 │ -118.584  │              │              │
╘════╧═════════════════╧═══════════╧═══════════╧══════════════╧══════════════╛

Draw final lens:

lens.draw(num_rays=5)