Lenses for telescope model?

[ppp-one]

Hi - loving the grand initiative here!

I'm trying to recreate a RC telescope, but have hit a road block with the lenses.

Is this possible to include in a SequentialSystem?

[byrdie]

Hi @ppp-one, thanks for giving my package a try and I really appreciate the feedback!

I've been focused on implementing reflecting telescopes so far, so I'm afraid that lenses aren't fully-supported yet, but they can be made to work.

My intent was for optika to model model lenses using the same procedure you would use for Zemax, where one surface represents the front of the lens and another represents the back surface.

I haven't quite fully thought through a material type for transmitting lenses since I haven't decided on a catalog for glass parameters. However, a very basic transmitting material (with a constant refractive index of 1.5 and no attenuation) can be represented using the following code:

import dataclasses

@dataclasses.dataclass(eq=False, repr=False)
class MyGlass(
    optika.materials.AbstractMaterial,
):
    shape = {}
    transformation = None
    is_mirror = False

    def index_refraction(
        self,
        rays: optika.rays.RayVectorArray,
    ) -> na.ScalarLike:
        return 1.5

    def attenuation(
        self,
        rays: optika.rays.RayVectorArray,
    ) -> na.ScalarLike:
        return 0 / u.mm

    def efficiency(
        self,
        rays: optika.rays.RayVectorArray,
        normal: na.AbstractCartesian3dVectorArray,
    ) -> na.ScalarLike:
        return 1

What optical constants are you using for your glass? This class can be extended as needed if you need to interpolate an array of refractive index measurements.

Through testing, I've discovered that my procedure for computing the field of view does not work if there are lenses between the pupil stop and the field stop. If this is the case for you, you will have to specify the field and manually trace rays through the system like in the following example:

aperture = optika.apertures.RectangularAperture(10 * u.mm)

front = optika.surfaces.Surface(
    name="front",
    sag=optika.sags.SphericalSag(100 * u.mm),
    material=MyGlass(),
    aperture=aperture,
    transformation=na.transformations.Cartesian3dTranslation(z=10 *u.mm),
    is_pupil_stop=True,
)
back = optika.surfaces.Surface(
    name="back",
    sag=optika.sags.SphericalSag(-100 * u.mm),
    aperture=aperture,
    transformation=na.transformations.Cartesian3dTranslation(z=15 * u.mm),
)

sensor = optika.surfaces.Surface(
    name="sensor",
    aperture=optika.apertures.RectangularAperture(1 * u.mm),
    is_field_stop=True,
    transformation=na.transformations.Cartesian3dTranslation(z=112.5 * u.mm),
)

pupil = na.Cartesian2dVectorLinearSpace(
    start=-1,
    stop=1,
    axis=na.Cartesian2dVectorArray("px", "py"),
    num=5,
    centers=True,
)

field = na.Cartesian2dVectorLinearSpace(
    start=-0.75 * u.deg,
    stop=0.75 * u.deg,
    axis=na.Cartesian2dVectorArray("fx", "fy"),
    num=3,
    centers=True,
)

wavelength = 500 * u.nm

grid_input = optika.vectors.ObjectVectorArray(
    wavelength=wavelength,
    field=field,
    pupil=pupil,
)

system = optika.systems.SequentialSystem(
    object=optika.surfaces.Surface(),
    surfaces=[
        front,
        back,
        sensor,
    ],
    grid_input=grid_input,
)

rays = system.raytrace(
    wavelength=wavelength,
    field=field,
    pupil=pupil,
    axis="surf",
    normalized_field=False,
)

with astropy.visualization.quantity_support():
    fig, ax = plt.subplots(constrained_layout=True)
    ax.set_aspect("equal")
    system.plot(
        ax=ax,
        components=("z", "y"),
        kwargs_rays=dict(
            color="tab:blue",
        ),
        color="black",
        zorder=10,
        plot_rays_vignetted=True,
        plot_rays=False,
    )
    na.plt.plot(
        rays.outputs.position,
        ax=ax,
        color="tab:blue",
        components=("z", "y"),
        axis="surf"
    )

Let me know how this works for you, I'm sure we can tweak it to suit your needs.

[ppp-one]

Thank you so much for the reply. I'm trying to recreate this set up:

Surface	Type	Comment	Radius	Thickness		Glass	Semi-Diameter		Conic
OBJ	Standard		Infinity	Infinity			Infinity		0.000000
1*	Standard		Infinity	1800.000000			509.887177		0.000000
*	Standard	M1	-4528.00000	-1665.38346	V	MIRROR	500.074962		-1.056959
3	Standard	M2	-1669.94912	1665.383456	P	MIRROR	142.369119		-3.786889
4	Standard		Infinity	330.390860	V		60.000000	U	0.000000
5*	Standard		335.200000	10.000019	V	F_SILICA	30.000000	U	0.000000
6*	Standard		409.730000	0.999995	V		30.000000	U	0.000000
7*	Standard		66.552000	8.027322	V	F_SILICA	30.000000	U	0.000000
8*	Standard		57.246000	100.581804	T		30.000000	U	0.000000
IMA	Standard		Infinity	-			21.954214		0.000000

This setup works on Zemax (however, my Zemax background is very limited at the moment). My end goal is see how the PSF changes for different focus position - figured it would be nice if I could control this programmatically directly in python than use zemax+the python api.

Here's my attempt, but I get ValueError: Max iterations exceeded. I suppose, maybe related to the issue you mentioned earlier? I set refractive index of the fused silica to 1.4585.

front = optika.surfaces.Surface(
    name="front",
    transformation=na.transformations.Cartesian3dTranslation(
        z = 0 * u.mm,
    ),
)

primary_mirror = optika.surfaces.Surface(
    name="primary_mirror",
    sag=optika.sags.ConicSag(
        radius=-4528 * u.mm,
        conic=-1.056959
    ),
    aperture=optika.apertures.CircularAperture(500.369119 * u.mm),
    material=optika.materials.Mirror(),
    is_pupil_stop=True,
    transformation=na.transformations.Cartesian3dTranslation(
        z = 1800 * u.mm,
    ),
)

secondary_mirror = optika.surfaces.Surface(
    name="secondary_mirror",
    sag=optika.sags.ConicSag(
        radius=-1669.94912 * u.mm,
        conic=-1.786889
    ),
    aperture=optika.apertures.CircularAperture(142.369119 * u.mm),
    material=optika.materials.Mirror(),
    transformation=na.transformations.Cartesian3dTranslation(
        z = (1800 - 1665.38346) * u.mm,
    ),
)

obscuration = optika.surfaces.Surface(
    name="obscuration",
    aperture=dataclasses.replace(secondary_mirror.aperture, inverted=True),
    transformation=secondary_mirror.transformation,
)

# lens 1
lens1front = optika.surfaces.Surface(
    name="lens1front",
    sag=optika.sags.SphericalSag(
        radius=335.2 * u.mm,
    ),
    aperture=optika.apertures.CircularAperture(30 * u.mm),
    material=MyGlass(),
    transformation=na.transformations.Cartesian3dTranslation(
        z = (1800 - 1665.38346 + 1665.38346 + 330.390860) * u.mm,
    ),
)

lens1back = optika.surfaces.Surface(
    name="lens1back",
    sag=optika.sags.SphericalSag(409.73 * u.mm),
    aperture=optika.apertures.CircularAperture(30 * u.mm),
    transformation=na.transformations.Cartesian3dTranslation(
        z=(1800 - 1665.38346 + 1665.38346 + 330.390860 + 10.000019 ) * u.mm),
)

# lens 2
lens2front = optika.surfaces.Surface(
    name="lens2front",
    sag=optika.sags.SphericalSag(
        radius=66.552 * u.mm,
    ),
    aperture=optika.apertures.CircularAperture(30 * u.mm),
    material=MyGlass(),
    transformation=na.transformations.Cartesian3dTranslation(
        z = (1800 - 1665.38346 + 1665.38346 + 330.390860 + 10.000019 + 0.99995) * u.mm,
    ),
)

lens2back = optika.surfaces.Surface(
    name="lens2back",
    sag=optika.sags.SphericalSag(
        radius= 57.246 * u.mm
    ),
    aperture=optika.apertures.CircularAperture(30 * u.mm),
    transformation=na.transformations.Cartesian3dTranslation(
        z=(1800 - 1665.38346 + 1665.38346 + 330.390860 + 10.000019 + 0.999995 + 8.027322) * u.mm
    ),
)


# define the imaging sensor surface
sensor = optika.sensors.ImagingSensor(
    name="sensor",
    width_pixel=10 * u.um,
    axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"),
    num_pixel=na.Cartesian2dVectorArray(1024, 1024),
    timedelta_exposure=1 * u.s,
    transformation=na.transformations.TransformationList([
        na.transformations.Cartesian3dTranslation(
            x=0 * u.mm,
            z=(1800 - 1665.38346 + 1665.38346 + 330.390860 + 10.000019 + 0.999995 + 8.027322 + 100.581804) * u.mm,
        ),
    ]),
    is_field_stop=True,
)

pupil = na.Cartesian2dVectorLinearSpace(
    start=-1,
    stop=1,
    axis=na.Cartesian2dVectorArray("px", "py"),
    num=5,
    centers=True,
)

field = na.Cartesian2dVectorLinearSpace(
    start=-0.1 * u.deg,
    stop=0.1 * u.deg,
    axis=na.Cartesian2dVectorArray("fx", "fy"),
    num=3,
    centers=True,
)

wavelength = 500 * u.nm

grid_input = optika.vectors.ObjectVectorArray(
    wavelength=wavelength,
    field=field,
    pupil=pupil,
)

system = optika.systems.SequentialSystem(
    object=optika.surfaces.Surface(),
    surfaces=[
        front,
        obscuration,
        primary_mirror,
        secondary_mirror,
        lens1front,
        lens1back,
        lens2front,
        lens2back,
    ],
    sensor=sensor,
    grid_input=grid_input,
)

rays = system.raytrace(
    wavelength=wavelength,
    field=field,
    pupil=pupil,
    axis="surf",
    normalized_field=False,
)

# plot the system
with astropy.visualization.quantity_support():
    fig, ax = plt.subplots(constrained_layout=True)
    ax.set_aspect("equal")
    system.plot(
        ax=ax,
        components=("z", "y"),
        kwargs_rays=dict(
            color="tab:blue",
        ),
        color="black",
        zorder=10,
        plot_rays_vignetted=True,
        plot_rays=False,
    )
    na.plt.plot(
        rays.outputs.position,
        ax=ax,
        color="tab:blue",
        components=("z", "y"),
        axis="surf"
    )

[byrdie]

I think this is related to the issue I mentioned earlier. To resolve this, I think it's best to also specify the pupil positions in physical coordinates for now

pupil = na.Cartesian2dVectorLinearSpace(
    start=-primary_mirror.aperture.radius,
    stop=primary_mirror.aperture.radius,
    axis=na.Cartesian2dVectorArray("px", "py"),
    num=5,
    centers=True,

and make sure to change normalized_pupil=False in the call to raytrace()

rays = system.raytrace(
    wavelength=wavelength,
    field=field,
    pupil=pupil,
    axis="surf",
    normalized_field=False,
    normalized_pupil=False,
)

[ppp-one]

Great 🕺, thank you so much. It rendered.

With this, I'm trying to follow the spot model tutorial here: https://optika.readthedocs.io/en/latest/tutorials/prime_focus.html

A copy and paste run yields this:

with the following error: ValueError: Max iterations exceeded again.

Changing angle = system.rayfunction_default.inputs.field.to(u.arcmin) to angle = rays.inputs.field.to(u.arcmin) removes this error, but yields a new error ValueError: the axisargument must beNoneor a subset of the broadcasted shape ofaandwhere, got {'fy'} for axis, but {'fx': 3} for shape

I don't think I fully understand the whole flow to fully debug.

[byrdie]

Agreed, I've had annoyances with the debugging process as well and I've been brainstorming ways to restructure things to make them easier to follow.

In this case, the error means that ScalarArray.mean() got called with an axis that doesn't exist since field is a sparse 2D grid. We can convert field into a dense grid using the broadcasted property, so your expression could be

angle = rays.inputs.field.broadcasted.to(u.arcmin)

but also, using what you know about the shape of your system, instead you could do

angle = rays.inputs.field.to(u.deg)
angle_x = angle.x
angle_y = angle.y

[byrdie]

Also, I noticed a small bug in this script where the vignetted rays are erroneously being included in the spot diagram. The call to scatter() should be

na.plt.scatter(
    position_relative.x,
    position_relative.y,
    ax=ax,
    s=5,
    where=rays.outputs.unvignetted,
)

[ppp-one]

That resolved that issue. The spot diagrams still look funny however.

I guess it's not constraining to the sensor. I tried my best to adjust the code, but without much success. In any case, this is what I tried:

I tried adapting the position_relative from the example, from:

position = system.rayfunction_default.outputs.position.to(u.um) # this yeilded an error

position_relative = position - position.mean(pupil.axes)

to (after trying to understand the source code of system.rayfunction_default):

axis = "_dummy"
_rays = rays.outputs

if sensor.transformation is not None:
    print("not none")
    _rays = sensor.transformation.inverse(_rays)

rays.outputs = _rays

position = rays.outputs.position.to(u.um)

position_relative = position - position.mean(pupil.axes)

However, this didn't help.

[byrdie]

Can you post the spot diagrams that Zemax outputs if you have it available? How small did you expect the spot size to be for this design?

I think there's either an issue with the optical prescription, or there's an unresolved issue with the internals of my refraction calculation, I will continue to investigate since this really should be working.

[ppp-one]

Using the same data as the original table (using C79-80 fused silica glass) produced this:

Hope it makes sense :). Thank you so much for your help.

[byrdie]

If it's not too much trouble, could you also post the same image but with the lenses disabled? I'm trying to isolate which part of the system is causing the issue, thanks!

[ppp-one]

Sure! With lenses removed (detector shifts)

Detector placed back to its place in the original design

[byrdie]

@ppp-one sorry for the delay, I think I've identified most of the problems.

The main reason that the spot diagrams looked wrong is that we were plotting the footprint for all surfaces, not the last surface, so we need to index the last element of the surf axis.

Another issue is that I think you have a typo, the conic constant of the secondary mirror should probably be something like -3.786889, not -1.786889.

With those fixes, the spot diagrams look much closer, although not exactly the same. I think the remaining discrepancy is due to a minor difference in the index of refraction. Let me know what you think, I'm pretty sure there's some way to extract the exact value of n Zemax is using to double check this hypothesis.

index = dict(surf=~0)
position = rays.outputs.position.to(u.um)[index]

position_relative = position - position.mean(pupil.axes)

with astropy.visualization.quantity_support():
    fig, ax = na.plt.subplots(
        axis_rows=field.axis.y,
        axis_cols=field.axis.x,
        nrows=field.num,
        ncols=field.num,
        sharex=True,
        sharey=True,
        figsize=(6, 6),
        constrained_layout=True,
    )
    na.plt.scatter(
        position_relative.x,
        position_relative.y,
        ax=ax,
        s=5,
        c=rays.outputs.unvignetted[index].astype(float),
    )

    ax_lower = ax[{field.axis.y: +0}]
    ax_upper = ax[{field.axis.y: ~0}]
    ax_left = ax[{field.axis.x: +0}]
    ax_right = ax[{field.axis.x: ~0}]

    na.plt.set_aspect("equal", ax=ax)
    na.plt.set_xlabel(f"$x$ ({position.x.unit:latex_inline})", ax=ax_lower)
    na.plt.set_ylabel(f"$y$ ({position.y.unit:latex_inline})", ax=ax_left)

    angle = rays.inputs.field.broadcasted.to(u.deg)
    angle_x = angle.x.mean(set(angle.axes) - {field.axis.x,})
    angle_y = angle.y.mean(set(angle.axes) - {field.axis.y,})
    na.plt.text(
        x=0.5,
        y=1,
        s=angle_x.to_string_array(),
        ax=ax_upper,
        transform=na.plt.transAxes(ax_upper),
        ha="center",
        va="bottom",
    )
    na.plt.text(
        x=1.05,
        y=0.5,
        s=angle_y.to_string_array(),
        ax=ax_right,
        transform=na.plt.transAxes(ax_right),
        ha="left",
        va="center",
    )
    fig.savefig("spot.png")

[ppp-one]

Amazing!! And good eye - really sorry about the typo.

There's one last thing I don't fully understand (perhaps I should open this as a separate issue - will do if you suggest):

• Why do some rays transmit through a surface if it's considered a mirror? e.g. through the secondary/primary?

I'm trying to recreate the spot diagram/PSF of donuts when the telescope is out of focus. The grand plan being - to help diagnose mirror alignment issues on a real telescope (I understand there are potential degeneracies here).

I have tried adding further obscurations (inverted and not), but I often get ValueError: stationary point detected. Any help beyond the amazing diagnosis/debugging you have done already is welcomed 🤞

[byrdie]

I'm glad we're getting closer! Thanks for being patient with this, the docs need to be much better.

Can you show me a picture of what you mean by rays transmitting through a surface?

Generally, the ValueError: stationary point detected is supposed to mean that it found at least one ray parallel to the surface. I'm still investigating this issue since it happens to me a lot, but so far only on optical systems which have an underlying problem, like a misplaced surface. In any case, I've been thinking that maybe it would be fine in this case to just mark these rays as vignetted and not halt the program.

[ppp-one]

Sure thing. For example, adding a new obscuration:

obscuration2 = optika.surfaces.Surface(
    name="obscuration2",
    aperture=dataclasses.replace(secondary_mirror.aperture),
    material=optika.materials.Mirror(),
    transformation=na.transformations.TransformationList(
        [
            # na.transformations.Cartesian3dRotationX((180) * u.deg),
            na.transformations.Cartesian3dTranslation(
                z=(1800 - 1665.38346 + 100) * u.mm,
                y=400.369119 * u.mm,
            ),
        ]
    ),
)

Rotating the mirror 180 degrees yeilds the same visual result.

Side note: Rotating the mirror to 10 degs does something weird:

[byrdie]

In the case of the obscuration, I would note that this is a sequential optical model and branching optical paths are not allowed, is that what you're trying to do here?

The reason the rays go through surfaces is that the code doesn't enforce that the ray intercepts must be positive, it is up to the user to make sure their optical system is actually sequential. I could be wrong, but doesn't Zemax work this way as well?

In the case of the second image, it is doing what I expect (which is the equal angle reflection law), it's just that the next surface (the primary) is behind the mirror which doesn't represent a physically-realizable sequential optical system.

[ppp-one]

I guess I don't fully understand the nature of the modelling here (or in Zemax for that matter).

I would just naively assume the ray tracing would follow from A to B, interacting with the materials it encounters. Neither of the scenarios above seem to make any physical sense to me.

[byrdie]

Ah, I think I may have misinterpreted what you were trying to do. This is supposed to be the surface representing the light blocked by the secondary mirror, correct? In that case you need to not set the material at all and set inverted to True.

My apologies, I thought this surface was a new thing you were introducing and I got confused.

[ppp-one]

Yeah, I've tried various combinations of invertion/placement/component order without success in attempting to achieve the (now) elusive donut shape (without the additional surface).

The last case of the 10 degree rotated mirror still very much confuses me - to me, it looks like it's behaving like a lens. Also, why is it then interacting with a primary mirror surface that also looks like it's rotated (not by me) -- I guess, a consequence of the sequential modelling?

[byrdie]

Yeah, I've tried various combinations of invertion/placement/component order without success in attempting to achieve the (now) elusive donut shape (without the additional surface).

I'm beginning to realize the donut shape you're describing is probably due to diffraction and cannot be modeled by raytracing. You likely need to use Zemax in Fourier optics mode, or a Python library like Poppy to model diffraction.

I would just naively assume the ray tracing would follow from A to B, interacting with the materials it encounters. Neither of the scenarios above seem to make any physical sense to me.

I think this is non-sequential raytracing which is exponentially harder than sequential which knows the order of the interacting surfaces in advance.

The last case of the 10 degree rotated mirror still very much confuses me - to me, it looks like it's behaving like a lens. Also, why is it then interacting with a primary mirror surface that also looks like it's rotated (not by me) -- I guess, a consequence of the sequential modelling?

This follows my expectation. Notice how the ray is bent at the new mirror by the equal-angle reflection law, it's just that the next surface in the sequence is to the right so the intercept found by optika is "inside" the mirror.

[ppp-one]

Thanks for directing me to Poppy. I'll have a deep dive.

As for the donut - I don't think it's diffraction caused (but I'm no optics expert), but rather it's the shadow of the secondary mirror that's visible. In the case above, when the secondary is shifted ~1mm (I think) from focus, you really see a donut.

e.g.