Feature Request: Add support for Gradient Index (GRIN) surfaces

[goldengrape]

Checklist

• I have searched the existing issues and discussions for a similar question or feature request.
• I have read the documentation and tried to find an answer there.
• I am using the latest version of Optiland.
• I have included all necessary context.

Thanks for taking the time to go through this — it really helps us help you!

Feature Request

Of course. Here is the final draft of the issue, written in English and incorporating all the details you've provided. You can copy and paste this directly into the GitHub issue submission form.

---

Feature Request

Is your feature request related to a problem? Please describe.
Currently, Optiland lacks support for surfaces with a gradient refractive index (GRIN), such as the gradient3 surface type available in Zemax. This limits the software's capability to design and simulate optical systems that utilize GRIN lenses. These lenses are crucial for correcting aberrations and reducing the number of optical elements in a system.

Furthermore, GRIN lenses are of significant importance in the fields of ophthalmology and optometry for simulating realistic eye models. Many established models, such as the Atchison and Liou-Brennan model eyes, rely on gradient refractive indices to accurately represent the human eye's crystalline lens. The absence of this feature in Optiland is a considerable limitation for researchers and engineers working on ophthalmic and vision science applications.

Describe the solution you'd like
I would like to request the addition of a new surface type that supports a gradient refractive index, similar to Zemax's gradient3. This would involve implementing a ray tracing algorithm capable of handling propagation through a medium with a continuously varying refractive index.

Describe alternatives you've considered
I have considered using a series of shells with uniform refractive indices to approximate a gradient index. However, this method is computationally intensive and would significantly slow down simulations, making it an impractical alternative for complex or iterative design processes.

Additional context
The importance of GRIN lenses is particularly evident in the simulation of human eye models, where they are used to accurately model the crystalline lens. For example, the Atchison (doi:10.1016/j.visres.2006.01.004) and Liou-Brennan model eyes (doi: 10.1364/josaa.14.001684) both utilize gradient-index lenses, which is a key reason for their anatomical and optical accuracy. Adding this functionality would make Optiland a much more powerful tool for ophthalmic and vision science research.

To assist with implementation, I am providing pseudocode for a ray tracing algorithm within a gradient medium. This algorithm uses the Runge-Kutta (RK4) method to solve the ray equation in a medium with a refractive index profile defined by radial and axial coefficients.

Here is the pseudocode:

// -------------------------------------------------------------------------
// Function: F_derivative(z, v, params)
// Purpose: Calculate the state derivative vector F = dv/dz at a given point (z, v)
// Inputs:
//   z: current z-coordinate (float)
//   v: current state vector [x, y, Tx, Ty, OP] (5D array)
//   params: coefficients of the Gradient3 medium (struct or object)
// Output:
//   Derivative vector F (5D array), or an error flag (e.g., null) if total internal reflection occurs
// -------------------------------------------------------------------------
FUNCTION F_derivative(z, v, params):
    // 1. Unpack the state vector
    x = v[0], y = v[1], Tx = v[2], Ty = v[3]

    // 2. Efficiently calculate n and its partial derivatives
    r_sq = x*x + y*y
    r_sq_2 = r_sq * r_sq
    r_sq_3 = r_sq_2 * r_sq

    n = params.n0 + params.nr2*r_sq + params.nr4*r_sq_2 + params.nr6*r_sq_3 +
        params.nz1*z + params.nz2*z*z + params.nz3*z*z*z

    g_r = 2*params.nr2 + 4*params.nr4*r_sq + 6*params.nr6*r_sq_2
    dn_dx = x * g_r
    dn_dy = y * g_r

    // 3. Calculate the Hamiltonian H and check for validity
    n_sq = n*n
    T_sq_transverse = Tx*Tx + Ty*Ty

    radicand = n_sq - T_sq_transverse // Value inside the square root
    IF radicand <= 0 THEN
        // Total internal reflection or ray is perpendicular to the z-axis, cannot continue tracing with z as the step
        RETURN null // Return an error flag
    END IF

    H = -sqrt(radicand)

    // 4. Calculate and return the derivative vector F
    inv_H = 1.0 / H // Pre-calculate the inverse to avoid multiple divisions
    F = [
        -Tx * inv_H,
        -Ty * inv_H,
        -n * dn_dx * inv_H,
        -n * dn_dy * inv_H,
        -n_sq * inv_H
    ]
    RETURN F
END FUNCTION


// -------------------------------------------------------------------------
// Main procedure: RayTrace_Gradient3_RK4
// -------------------------------------------------------------------------
PROCEDURE RayTrace_Gradient3_RK4:
    // 1. Initialization
    // Define medium parameters
    gradient_params = {n0, nr2, nr4, nr6, nz1, nz2, nz3}

    // Define initial ray state
    z_start = 0.0
    z_end = 10.0
    step_size = 0.1

    initial_pos = [x0, y0]
    initial_dir = [dir_x0, dir_y0, dir_z0] // Unit direction vector

    // Calculate the initial state vector based on initial conditions
    n_start = calculate_n(z_start, initial_pos, gradient_params) // Requires a helper function
    Tx0 = n_start * initial_dir[0]
    Ty0 = n_start * initial_dir[1]

    z = z_start
    v = [initial_pos[0], initial_pos[1], Tx0, Ty0, 0.0] // [x, y, Tx, Ty, OP]

    // Store the trajectory
    ray_path = []
    ADD {z: z, state: v} TO ray_path

    // 2. RK4 main loop
    WHILE z < z_end:
        // Ensure the last step ends exactly at z_end
        current_step = min(step_size, z_end - z)

        // Calculate K1
        K1 = F_derivative(z, v, gradient_params)
        IF K1 IS null THEN BREAK // Error handling

        // Calculate K2
        v_temp2 = v + (current_step/2.0) * K1
        K2 = F_derivative(z + current_step/2.0, v_temp2, gradient_params)
        IF K2 IS null THEN BREAK // Error handling

        // Calculate K3
        v_temp3 = v + (current_step/2.0) * K2
        K3 = F_derivative(z + current_step/2.0, v_temp3, gradient_params)
        IF K3 IS null THEN BREAK // Error handling

        // Calculate K4
        v_temp4 = v + current_step * K3
        K4 = F_derivative(z + current_step, v_temp4, gradient_params)
        IF K4 IS null THEN BREAK // Error handling

        // 3. Update the state vector and z-coordinate
        v = v + (current_step / 6.0) * (K1 + 2*K2 + 2*K3 + K4)
        z = z + current_step

        ADD {z: z, state: v} TO ray_path

    END WHILE

    // 4. Output the results
    IF z >= z_end THEN
      PRINT "Ray tracing completed successfully."
      PRINT "Final state at z = ", z, ":"
      PRINT "  Position (x, y) = (", v[0], ", ", v[1], ")"
      PRINT "  Direction (Tx, Ty) = (", v[2], ", ", v[3], ")"
      PRINT "  Optical Path (OP) = ", v[4]
    ELSE
      PRINT "Ray tracing terminated early due to an error (e.g., total internal reflection) at z = ", z
    END IF

END PROCEDURE

[HarrisonKramer]

Hi @goldengrape,

Thanks for opening this feature request and for sharing the detailed background and pseudocode. Supporting GRIN media would indeed open up valuable applications, and it's a feature we've wanted to add for a while. Right now, Optiland’s tracing model is surface-to-surface: rays travel in straight lines between boundaries. GRIN propagation will of course require a different approach, where the trajectory is integrated step by step through a continuously varying medium. This makes it more challenging to implement, since it involves numerical integration and new handling of boundaries. I need to investigate how best this might fit into the current framework. I am open to suggestions if you have any.

Your pseudocode gives a clear starting point for how such an algorithm could look, and it does fit conceptually into Optiland’s workflow. We’ll add GRIN support to the roadmap, though I don't have a timeline yet for when this will be implemented. If anyone is interested in contributing, this would be a great project to collaborate on, and we’d be glad to provide guidance.

I will keep this issue open to track the progress of this work. Thanks again!

Kramer

[goldengrape]

I attempted to implement a GRIN lens using a shell-based discretization approach. Beyond the inherent complexity of the programmatic design, a critical limitation emerged under sequential ray tracing: rays cannot traverse the peripheral regions of a small-aperture lens to reach the next surface. This constraint effectively prevents ray propagation within the shell structure. Consequently, the alternative strategy of approximating the GRIN medium with concentric shells of uniform refractive index proves infeasible. For further context, refer to the discussion on #345

[HarrisonKramer]

Hi @goldengrape,

Thanks for working on a proposal. I looked through your code and I recognize the difficulties.

While I have not yet checked how best to implement this, I have a few thoughts on a strategy.

• We could implement a new material type, e.g., GRINMaterial, which subclasses BaseMaterial
• Currently, rays propagate through a medium in RealRays - this is just straight line propagation.
• We can intervene here, instead passing the ray propagation to a new method in GRINMaterial. We could use your numerical integration approach here to identify the path through the material.
• The propagation continues until the ray(s) exit the surface boundary.
• There are other subtle details, like how to track/save path through the material, especially for later analysis or visualization, but these aren't critical to handle just yet.

While I think something like this would work, it would require a relatively large refactor. It may also lead to overly complex code that is less maintainable, so we'd need to be careful in how we decide to implement it.

I will think further about this, as it's a feature I'd love to have in the package. If you have further thoughts, please feel free to share.

Best,
Kramer

[goldengrape]

I've forked the repository (https://github.com/goldengrape/optiland) and am currently working on implementing GRIN surface support. Broadly speaking, I need to add components across three modules to enable GRIN functionality: Geometry (Surface), Physics (Material), and Behavior (Propagation). I may also need to modify optic/optic.py. ( Plan: https://github.com/goldengrape/optiland/blob/master/docs/GRIN_design/GRIN_design_and_implementation.rst )

The work is still in progress, and I plan to submit a pull request once I’ve completed testing. My programming skills are limited, so this task is quite challenging for me, hopefully I’ll be able to complete it. I'm relying heavily on AI assistance (Gemini), that's why I use 'icontract', but I will remove them before pull request.

[HarrisonKramer]

Hi @goldengrape,

Thanks for sharing your progress and detailed plan! This will be a solid addition to the codebase.

Before we dive into the full implementation, I think we should take a step back and look at how this fits into the overall structure. Right now, Optiland’s tracing logic assumes straight-line propagation between surfaces. GRIN tracing changes that completely, so we’ll need to make sure the core architecture can handle both without turning into a tangle of special cases. Of course, it would be possible to already implement the GRIN propagation functionality, but it will almost certainly lead to a poor code design that will be hard to maintain and extend.

In your plan, the tracer checking for something like isinstance(surface, GradientBoundarySurface) would technically work, but it ties the core logic directly to one specific feature. That’s risky long-term. Instead, I think we should introduce a small abstraction layer - something like a “Propagation Model.”

We could move the current straight-line propagation into e.g., a HomogeneousPropagation class, and your RK4-based GRIN solver could live in a GrinPropagation class. The main tracer would then just call propagate() on whichever model applies, with no hardcoded checks. That keeps things modular and makes it easy to add new propagation types later, as all propagation is abstracted away.

So, I’d suggest we focus first on that refactor. Once the propagation interface exists, your GRIN code can plug right in with minimal friction.

You’ve done excellent groundwork here, and it’s already helping shape the direction of the project. Let’s use this opportunity to get the foundation right so GRIN (and future models) can fit in cleanly. I will take a day or two to think through this new "propagation" model and try to implement something. I'm confident it can be done with minimal friction and that will make it much easier for you to be able to impement the GRIN functionality.

Thanks again!
Kramer

[goldengrape]

I have updated the plan according to your arrangement:
https://github.com/goldengrape/optiland/blob/master/docs/GRIN_design/GRIN_design_and_implementation.rst
https://github.com/goldengrape/optiland/blob/master/docs/GRIN_design/todolist.md

I have separately implemented interaction/gradient_propagation, materials/gradient_material, and surfaces/gradient_surface.

I plan to submit a pull request after completing phase 2 in the todolist, because phase 3 might require modifications to optic/optic.

[HarrisonKramer]

Hey @goldengrape,

Thanks for the detailed proposal! I've spent quite a bit of time reviewing your docs and thinking this through in more detail. Most of your proposal is solid, but I would suggest some slight variations. Let me know what you think.

I think the real issue is the current PropagationModel API. Right now it's propagate(rays, t), which assumes the surface pre-calculates a straight-line distance. This is exactly what breaks for GRIN.

Here's what I'm thinking instead:

1. Change the API signature from propagate(rays, t) to e.g., propagate(rays, exit_geometry). This lets the propagation model control the entire journey, not just execute a pre-calculated step.

2. Simplify Surface.trace - remove the distance calculation and just call self.material_pre.propagation_model.propagate(rays, self.geometry). The surface shouldn't care how rays arrive, just handle the interaction afterward.

3. Update HomogeneousPropagation to implement the new interface (it'll do the distance calc internally, same logic as before).

4. Create GrinPropagation - new class that does the RK4 integration, calling e.g., material.get_index_and_gradient() in a loop and checking exit_geometry to know when to stop.

5. Create GradientMaterial - instantiates its own GrinPropagation in __init__ and provides the get_index_and_gradient() method.

6. Visualization Updates - as mentioned previously, we will also need to change how visualization is performed. This is not critical for now, but we will need to update this. I have an idea on how to tackle this, e.g., by using an optional recorder argument that records the (x, y, z) positions during the raytrace when it's passed, but that can come later.

One big benefit is that we don't touch SurfaceGroup.trace at all, and adding new propagation types in the future is just a matter of creating new Material/Propagation pairs. We also don't need a new GRIN Surface type. Much cleaner for long-term maintainability.

Just some thoughts so far, as I'm not 100% certain on all details. Let me know what you think, or if you disagree on anything. I am happy to iterate. I think your proposal and your code is quite consistent with this approach.

Regards,
Kramer