Giter VIP home page Giter VIP logo

Comments (4)

kklmn avatar kklmn commented on August 17, 2024

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.

from xrt.

HGL2023 avatar HGL2023 commented on August 17, 2024

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.
picture

from xrt.

HGL2023 avatar HGL2023 commented on August 17, 2024

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!
fig 1
fig 2

from xrt.

HGL2023 avatar HGL2023 commented on August 17, 2024

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

from xrt.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.