Some doubts about custom faceted lenses

[HGL2023]

Hi, there are a couple of questions about lenses with custom face shapes that I have not been able to figure out, so I'd appreciate your guidance.
1, I need to set my lens shape to "DoubleParaboloidLens", Do I need to rewrite cl_local_z and cl_local_n as custom face functions? as shown in the figure below .
(I think the XRT program relies on the text content in cl_local_z and cl_local_n to determine the lens shape, but I think the XRT program relies on the text content in cl_local_z and cl_local_n to determine the lens shape. cl_local_n to determine the lens shape.

2, I follow the tutorial to add a new module, in the new module to add new variables and change the face type function equation, successfully imported into xrt.qook, the generated script runs perfectly, but the program will prompt "Incorrect parameters , class customLens01 is not initialized" as shown in Figure 2 below. Is there a good solution for this situation.


Thank you for your guidance！

[kklmn]

Hi,

You are lumping together several problems.

First, don't use OpenCL (targetOpenCL=None, as by default). Then you don't need cl_ prefixed functions. You need to modify local_z1() and local_n1() for the upstream lens surface and local_z2() and local_n2() for the downstream one. You may also want to modify get_nCRL() to get a correct number of lenses for a given focal length, which depends on whether your lens focalizes in both directions and by both surfaces.

You can try your lens class with one of the examples in examples/withRaycing/04_Lenses.

If you want to boost the speed with a decent GPU, you specify targetOpenCL as needed. Then you also need your custom implementation of those OpenCL functions.

When the class is functional, whether with or without OpenCL connection, it should be importable into xrtQook.

[HGL2023]

Hello.

In XRT originally z = (x2 + y2) / (4 * self.focus), but I changed the formula for z. The original opening of my formula was downwards, and the program kept getting errors after replacing it directly, but adding a negative sign in front of the formula (which puts the opening upwards to keep the same direction as the opening of the original formula in XRT) made it work fine.

In this case, I think I need to make some changes in local_z2 to realize my lens correctly, but I don't really know how to do it, so I'd appreciate your guidance!

The red box in the picture is where I added the minus sign in front of my original equation.

[HGL2023]

Hello!
Thanks a lot for the previous instructions!
Recently, I've come across a new problem, our lenses are defined by q1,q2 for the front side and q3,q4 for the back side, and the number of lenses is input by myself, which means it's not simply the same lens for both front and back sides (the lenses are shown in Figure 1 below). So I defined local_z1, local_n1 and local_z2, local_n2 (I uploaded them) but the program eventually ran with an error (Figure 2), which left me at a bit of a loss as to how to modify it. (I think there are some steps I missed while modifying, but I just can't find them)
I hope I can get your help, thanks!

[HGL2023]

from xrt.backends.raycing.oes import *

class CustomLens(Plate):
hiddenMethods = Plate.hiddenMethods + ['double_refract']
cl_plist = ("zmax", "focus", "q1", "q2")
cl_local_z = """
float local_z(float8 cl_plist, int i, float x, float y)
{
float res;
res = -(-887444.2084046695 * (((x ** 2 + y ** 2 + 36594846827921.74) ** 0.5) - 6049367.47337453));
if (res > cl_plist.s0) res = cl_plist.s0;
return res;
}"""
cl_local_n = """
float3 local_n(float8 cl_plist, int i, float x, float y)
{
float3 res;
res.s0 = -(887444.2084046695 * (x ** 2 + y ** 2 + 36594846827921.74) ** (-0.5) * x);
res.s1 = -(887444.2084046695 * (x ** 2 + y ** 2 + 36594846827921.74) ** (-0.5) * y);
float z = -(-887444.2084046695 * (((x ** 2 + y ** 2 + 36594846827921.74) ** 0.5) - 6049367.47337453));
if (z > cl_plist.s0)
{
res.s0 = 0;
res.s1 = 0;
}
res.s2 = 1.;
return normalize(res);
}"""

def __init__(self, *args, **kwargs):
    kwargs = self.__pop_kwargs(**kwargs)
    Plate.__init__(self, *args, **kwargs)

@property
def nCRL(self):
    return self._nCRL

@nCRL.setter
def nCRL(self, nCRL):
    if isinstance(nCRL, (int, float)):
        self._nCRL = max(int(round(nCRL)), 1)
        self._nCRLlist = None
    elif isinstance(nCRL, (list, tuple)):
        self._nCRL = max(int(round(self.get_nCRL(*nCRL))), 1)
        self._nCRLlist = copy.copy(nCRL)
    #            print('nCRL={0}'.format(nCRL))
    else:
        self._nCRL = 1

#            raise ValueError("wrong nCRL value!")
@property                                                              ///four new variables q1,q2,q3,q4 are added here///
def q1(self):
    return self._q1
@q1.setter
def q1(self, q1):
    self._q1 = q1
    if hasattr(self, '_nCRLlist'):
        if self._nCRLlist is not None:
            self.nCRL = self._nCRLlist

@property
def q2(self):
    return self._q2
@q2.setter
def q2(self, q2):
    self._q2 = q2
    if hasattr(self, '_nCRLlist'):
        if self._nCRLlist is not None:
            self.nCRL = self._nCRLlist

@property
def q3(self):
    return self._q3

@q3.setter
def q3(self, q3):
    self._q3 = q3
    if hasattr(self, '_nCRLlist'):
        if self._nCRLlist is not None:
            self.nCRL = self._nCRLlist

@property
def q4(self):
    return self._q4

@q4.setter
def q4(self, q4):
    self._q4 = q4
    if hasattr(self, '_nCRLlist'):
        if self._nCRLlist is not None:
            self.nCRL = self._nCRLlist
@property
def focus(self):
    return self._focus
@focus.setter
def focus(self, focus):
    self._focus = focus
    if hasattr(self, '_nCRLlist'):
        if self._nCRLlist is not None:
            self.nCRL = self._nCRLlist

def __pop_kwargs(self, **kwargs):
    self.focus = kwargs.pop('focus', None)
    self.zmax = kwargs.pop('zmax', None)
    self.nCRL = kwargs.pop('nCRL', 1)
    self.q1 = kwargs.pop('q1', None)
    self.q2 = kwargs.pop('q2', None)
    self.q3 = kwargs.pop('q3', None)
    self.q4 = kwargs.pop('q4', None)
    kwargs['pitch'] = kwargs.get('pitch', np.pi / 2)
    return kwargs

def assign_auto_material_kind(self, material):
    material.kind = 'lens'

def local_z1(self, x, y):
    """Determines the normal vector of OE at (x, y) position."""
    # print(self.q1, self.q2)
    if self.q1 is None:
        z = -(-887444.2084046695 * (((x ** 2 + y ** 2 + 36594846827921.74) ** 0.5) - 6049367.47337453))
    else:
        z = -(-2 * (((self.q1 ** 2 - self.q1 * self.q2 + self.q2 ** 2) / 9) ** 0.5) * np.cos(np.arccos((((self.q1 + self.q2) * (9 * self.q1 * self.q2 - 2 * (self.q1 + self.q2) ** 2) / 54 - (self.q1 - self.q2) * (x ** 2 + y ** 2) / 4 * (-1.12683267e-6)) / ((self.q1 ** 2 - self.q1 * self.q2 + self.q2 ** 2) / 9) ** 1.5) / 3)) + (self.q1 + self.q2) / 3)
    if self.zmax is not None:
        z[z > self.zmax] = self.zmax
    return z

def local_z2(self, x, y):
    # if self.q3 is not None:
    z = -(-2 * (((self.q3 ** 2 - self.q3 * self.q4 + self.q4 ** 2) / 9) ** 0.5) * np.cos(np.arccos((((self.q3 + self.q4) * (9 * self.q3 * self.q4 - 2 * (self.q3 + self.q4) ** 2) / 54 - (self.q3 - self.q4) * (x ** 2 + y ** 2) / 4 * (-1.12683267e-6)) / ((self.q3 ** 2 - self.q3 * self.q4 + self.q4 ** 2) / 9) ** 1.5) / 3)) + (self.q3 + self.q4) / 3)
    if self.zmax is not None:
        z[z > self.zmax] = self.zmax
    return z

def local_n1(self, x, y):
    # just flat:
    if self.q1 is None:
        a = -0.1467003306231*x/(2.73262518272656e-14*x**2 + 2.73262518272656e-14*y**2 + 1)**0.5  # -dz/dx   a = 887444.2084046695 * (x ** 2 + y ** 2 + 36594846827921.74) ** (-0.5) * x
        b = -0.1467003306231*y/(2.73262518272656e-14*x**2 + 2.73262518272656e-14*y**2 + 1)**0.5  # -dz/dy   b = 887444.2084046695 * (x ** 2 + y ** 2 + 36594846827921.74) ** (-0.5) * y
    else:
        a = -3.7561089e-7 * x * (self.q1 - self.q2) * np.sin(0.333 * np.arccos((2.817081675e-7 * (self.q1 - self.q2) * (x ** 2 + y ** 2) + (self.q1 + self.q2) * (9 * self.q1 * self.q2 - 2 * (self.q1 + self.q2) ** 2) / 54) / (self.q1 ** 2 / 9 - self.q1 * self.q2 / 9 + self.q2 ** 2 / 9) ** 1.5)) / (np.sqrt(-(1.5212241045e-5 * (self.q1 - self.q2) * (x ** 2 + y ** 2) + (self.q1 + self.q2) * (9 * self.q1 * self.q2 - 2 * (self.q1 + self.q2) ** 2)) ** 2 / (2916 * (self.q1 ** 2 / 9 - self.q1 * self.q2 / 9 + self.q2 ** 2 / 9) ** 3.0) + 1) * (self.q1 ** 2 / 9 - self.q1 * self.q2 / 9 + self.q2 ** 2 / 9) ** 1.0)
        b = -3.7561089e-7 * y * (self.q1 - self.q2) * np.sin(0.333 * np.arccos((2.817081675e-7 * (self.q1 - self.q2) * (x ** 2 + y ** 2) + (self.q1 + self.q2) * (9 * self.q1 * self.q2 - 2 * (self.q1 + self.q2) ** 2) / 54) / (self.q1 ** 2 / 9 - self.q1 * self.q2 / 9 + self.q2 ** 2 / 9) ** 1.5)) / (np.sqrt(-(1.5212241045e-5 * (self.q1 - self.q2) * (x ** 2 + y ** 2) + (self.q1 + self.q2) * (9 * self.q1 * self.q2 - 2 * (self.q1 + self.q2) ** 2)) ** 2 / (2916 * (self.q1 ** 2 / 9 - self.q1 * self.q2 / 9 + self.q2 ** 2 / 9) ** 3.0) + 1) * (self.q1 ** 2 / 9 - self.q1 * self.q2 / 9 + self.q2 ** 2 / 9) ** 1.0)

    if self.zmax is not None:
        if self.q1 is None:
            z = -(-887444.2084046695 * ((x ** 2 + y ** 2 + 36594846827921.74) ** 0.5 - 6049367.47337453))
        else:
            z = -(-2 * (((self.q1 ** 2 - self.q1 * self.q2 + self.q2 ** 2) / 9) ** 0.5) * np.cos(np.arccos((((self.q1 + self.q2) * (9 * self.q1 * self.q2 - 2 * (self.q1 + self.q2) ** 2) / 54 - (self.q1 - self.q2) * (x ** 2 + y ** 2) / 4 * (-1.12683267e-6)) / ((self.q1 ** 2 - self.q1 * self.q2 + self.q2 ** 2) / 9) ** 1.5) / 3)) + (self.q1 + self.q2) / 3)
        if isinstance(a, np.ndarray):
            a[z > self.zmax] = 0
        if isinstance(b, np.ndarray):
            b[z > self.zmax] = 0
    c = np.ones_like(x)
    norm = (a ** 2 + b ** 2 + 1) ** 0.5
    return [a / norm, b / norm, c / norm]

def local_n2(self, x, y):
    if self.q3 is not None:
        a = -3.7561089e-7 * x * (self.q3 - self.q4) * np.sin(0.333 * np.arccos((2.817081675e-7 * (self.q3 - self.q4) * (x ** 2 + y ** 2) + (self.q3 + self.q4) * (9 * self.q3 * self.q4 - 2 * (self.q3 + self.q4) ** 2) / 54) / (self.q3 ** 2 / 9 - self.q3 * self.q4 / 9 + self.q4 ** 2 / 9) ** 1.5)) / (np.sqrt(-(1.5212241045e-5 * (self.q3 - self.q4) * (x ** 2 + y ** 2) + (self.q3 + self.q4) * (9 * self.q3 * self.q4 - 2 * (self.q3 + self.q4) ** 2)) ** 2 / (2916 * (self.q3 ** 2 / 9 - self.q3 * self.q4 / 9 + self.q4 ** 2 / 9) ** 3.0) + 1) * (self.q3 ** 2 / 9 - self.q3 * self.q4 / 9 + self.q4 ** 2 / 9) ** 1.0)
        b = -3.7561089e-7 * y * (self.q3 - self.q4) * np.sin(0.333 * np.arccos((2.817081675e-7 * (self.q3 - self.q4) * (x ** 2 + y ** 2) + (self.q3 + self.q4) * (9 * self.q3 * self.q4 - 2 * (self.q3 + self.q4) ** 2) / 54) / (self.q3 ** 2 / 9 - self.q3 * self.q4 / 9 + self.q4 ** 2 / 9) ** 1.5)) / (np.sqrt(-(1.5212241045e-5 * (self.q3 - self.q4) * (x ** 2 + y ** 2) + (self.q3 + self.q4) * (9 * self.q3 * self.q4 - 2 * (self.q3 + self.q4) ** 2)) ** 2 / (2916 * (self.q3 ** 2 / 9 - self.q3 * self.q4 / 9 + self.q4 ** 2 / 9) ** 3.0) + 1) * (self.q3 ** 2 / 9 - self.q3 * self.q4 / 9 + self.q4 ** 2 / 9) ** 1.0)
    if self.zmax is not None:
        z = -(-2 * (((self.q3 ** 2 - self.q3 * self.q4 + self.q4 ** 2) / 9) ** 0.5) * np.cos(np.arccos((((self.q3 + self.q4) * (9 * self.q3 * self.q4 - 2 * (self.q3 + self.q4) ** 2) / 54 - (self.q3 - self.q4) * (x ** 2 + y ** 2) / 4 * (-1.12683267e-6)) / ((self.q3 ** 2 - self.q3 * self.q4 + self.q4 ** 2) / 9) ** 1.5) / 3)) + (self.q3 + self.q4) / 3)

        if isinstance(a, np.ndarray):
            a[z > self.zmax] = 0
        if isinstance(b, np.ndarray):
            b[z > self.zmax] = 0
    c = np.ones_like(x)
    norm = (a ** 2 + b ** 2 + 1) ** 0.5
    return [a / norm, b / norm, c / norm]
                                                              ////the following program has not been modified////
def get_nCRL(self, f, E):
    nCRL = 1
    if all([hasattr(self, val) for val in ['focus', 'material']]):
        if self.focus is not None and self.material is not None:
            if isinstance(self, DoubleParaboloidLens):
                nFactor = 0.5
            elif isinstance(self, ParabolicCylinderFlatLens):
                nFactor = 2.
            else:
                nFactor = 1.
            nCRL = 2 * self.focus / float(f) / \
                   (1. - self.material.get_refractive_index(E).real) * nFactor
    return nCRL

def multiple_refract(self, beam=None, needLocal=True,
                     returnLocalAbsorbed=None):
    """
    Sequentially applies the :meth:`double_refract` method to the stack of
    lenses, center of each of *nCRL* lens is shifted by *zmax* mm
    relative to the previous one along the beam propagation direction.
    Returned global beam emerges from the exit surface of the last lens,
    returned local beams correspond to the entrance and exit surfaces of
    the first lens.

    *returnLocalAbsorbed*: None, 0 or 1
        If not None, returns the absorbed intensity in local beam. If
        equals zero, the last local beam returns total absorbed intensity,
        otherwise the absorbed intensity on single element of the stack.


    .. Returned values: beamGlobal, beamLocal1, beamLocal2
    """
    if self.bl is not None:
        self.bl.auto_align(self, beam)
    if self.nCRL == 1:
        self.centerShift = np.zeros(3)
        return self.double_refract(beam, needLocal=needLocal,
                                   returnLocalAbsorbed=returnLocalAbsorbed)
    else:
        tmpFlowSource = self.bl.flowSource
        self.bl.flowSource = 'multiple_refract'
        tempCenter = [c for c in self.center]
        beamIn = beam
        step = 2. * self.zmax + self.t \
            if isinstance(self, DoubleParaboloidLens) else self.zmax + self.t
        for ilens in range(self.nCRL):
            if isinstance(self, ParabolicCylinderFlatLens):
                self.roll = -np.pi / 4 if ilens % 2 == 0 else np.pi / 4
            lglobal, tlocal1, tlocal2 = self.double_refract(
                beamIn, needLocal=needLocal)
            if self.zmax is not None:
                toward = raycing.rotate_point(
                    [0, 0, 1], self.rotationSequence, self.pitch,
                    self.roll + self.positionRoll, self.yaw)
                self.center[0] -= step * toward[0]
                self.center[1] -= step * toward[1]
                self.center[2] -= step * toward[2]
            beamIn = lglobal
            if ilens == 0:
                llocal1, llocal2 = tlocal1, tlocal2
        self.centerShift = step * np.array(toward)
        self.center = tempCenter
        self.bl.flowSource = tmpFlowSource

        if returnLocalAbsorbed is not None:
            if returnLocalAbsorbed == 0:
                absorbedLb = rs.Beam(copyFrom=llocal2)
                absorbedLb.absorb_intensity(beam)
                llocal2 = absorbedLb
            elif returnLocalAbsorbed == 1:
                absorbedLb = rs.Beam(copyFrom=llocal1)
                absorbedLb.absorb_intensity(beam)
                llocal1 = absorbedLb
            elif returnLocalAbsorbed == 2:
                absorbedLb = rs.Beam(copyFrom=llocal2)
                absorbedLb.absorb_intensity(llocal1)
                llocal2 = absorbedLb
        llocal2.parent = self
        raycing.append_to_flow(self.multiple_refract,
                                [lglobal, llocal1, llocal2],
                                inspect.currentframe())
        return lglobal, llocal1, llocal2