Ray traces hit fictitious surfaces when using aspheres

[ElliotB256]

Hello,

I'm attempting to use raypier to simulate some coupling optics based on aspheres. I had some surprising results, and after visualising I can see that the rays seem to ignore the physical extent of the asphere. In the image below, you can see the rays refracting off non-existent surfaces:

I have defined my aspheres using a definition analogous to those used in the aspheres module for other Thorlabs lenses (see, for example, the definition here):

class TL_354330(AsphericLens):
    name = "354330"
    diameter = 5
    CT = 3.214
    A_curvature = 1/-3.669210850149012300E-001
    A_conic = 5.421234096046E-1 * -1
    A4 = 2.077710281508E-5* 1
    A6 = 1.373929859525E-4* 1
    A8 = -2.029327872896E-5* 1
    A10 = 2.961194418229E-6* 1
    A12 = 0
    A14 = 0
    A16 = 0
    B_curvature = 1/3.356784405862046600E-001
    B_conic = 1.210285287639E+1* -1
    B4 = -8.068904748567E-3* 1
    B6 = 1.725894132603E-3* 1
    B8 = -1.839071992368E-4* 1
    B10 = 8.936267158057E-6* 1
    B12 = 0
    B14 = 0
    B16 = 0

The phantom refractions I am seeing clearly trace the continued surface of the asphere past the diameter - it is as if the diameter is only used for drawing the mesh and ignored for ray tracing:

[ElliotB256]

The behavior is also visible in the check_thorlabs_aspheric example configuration. Here I have modified the ray source to have only one ray per ring, and a larger number of rings, so I can see a fan of rays from the side:

Interaction of rays outside the extent of the optic is visible:

This can lead to surprising effects far from the optics, presumably as a result of the asphere function hitting some divergence:

[ElliotB256]

I'm also seeing the behavior on PlanoConicLens but not for PlanoConvexLens lenses, which leads me to assume the issue is in SurfaceOfRotationLens

[ElliotB256]

CircularFace seems to explicitly test diameter in core/cfaces.pyx:

        X = p1.x + h*(p2.x-p1.x) - self.offset
        Y = p1.y + h*(p2.y-p1.y)
        if is_base_ray and (X*X + Y*Y) > (d*d/4):
            #print "X", X, "Y", Y
            return NO_INTERSECTION

ConicRevolutionFace uses the following routine:

        if is_base_ray and not (<Shape>(self.shape)).point_inside_c(pt1.x, pt1.y):
            return NO_INTERSECTION

As does AsphericFace:

        if is_base_ray and not (<Shape>(self.shape)).point_inside_c(pt1.x, pt1.y):
            return NO_INTERSECTION

Both of these faces inherit from ShapedFace, which specifies the shape parameter:

    def __cinit__(self, **kwds):
        self.shape = kwds.get("shape", Shape())

However, as far as I can see nowhere in the Aspheric or PlanoConicLens code ever actually sets the shape. For example, this is how faces are defined in AsphericLens:

AsphericFace(owner=self, diameter=self.diameter,
                                material=self.material,
                                conic_const=self.A_conic,
                                z_height=0.0, curvature=self.A_curvature,
                                A4=-self.A4, A6=-self.A6, A8=-self.A8, A10=-self.A10,
                                A12=-self.A12, A14=-self.A14, A16=-self.A16,
                                invert_normals=True,
                                atol=1e-16)

My understanding is that no shape is set for these lenses, and so they have no extent.

[ElliotB256]

After making the following change to AsphericFace's definition so that it includes a circular shape, I find that the ray tracing works as expected:

    def make_faces(self):
        fl = [AsphericFace(owner=self, diameter=self.diameter,
                                material=self.material,
                                conic_const=self.A_conic,
                                z_height=0.0, curvature=self.A_curvature,
                                A4=-self.A4, A6=-self.A6, A8=-self.A8, A10=-self.A10,
                                A12=-self.A12, A14=-self.A14, A16=-self.A16,
                                invert_normals=True,
                                shape=cCircleShape(radius=self.diameter/2),
                                atol=1e-16), 
                AsphericFace(owner=self, diameter=self.diameter,
                                material=self.material,
                                conic_const=self.B_conic,
                                z_height=self.CT, curvature=self.B_curvature,
                                A4=-self.B4, A6=-self.B6, A8=-self.B8, A10=-self.B10,
                                A12=-self.B12, A14=-self.B14, A16=-self.B16,
                                shape=cCircleShape(radius=self.diameter/2),
                                atol=1e-16)]

[ElliotB256]

(If you are happy with this fix, I can raise a PR to fix AsphericFace and ConicRevolutionFace. It's probably sensible to put this into SurfaceOfRotationLens).

[bryancole]

I guess my email replies didn't make it to github so we ended up fixing this in parallel. Anyway, I've pushed some changes which are equivalent to yours (on Master). Hence, closing this issue. Thanks for the input and finding the bug. BTW. SurfaceOfRotationLens is an abstract base-class so it doesn't define any faces, hence there's nowhere to put a Shape object in it.

[ElliotB256]

Thanks for the resolution!