Comments (4)
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.
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.
from xrt.
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!
from xrt.
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)
- simulation of VLS HOT 6
- build spherical mirror from elliptical optic HOT 4
- simulation of Hyperboloid mirror HOT 3
- Question about plotting energy spectrum HOT 6
- How to change the screen size or resolution when viewing a small focus image? HOT 2
- problem with xrtBentXtal in Laue geometry HOT 5
- plot error with matplotlib 3.8.2 HOT 6
- ParabolicalMirrorParam optics with figure error HOT 3
- xyz_to_param and param_to_xyz for class ParabolicalMirrorParam HOT 2
- How to simulate nested shells HOT 9
- Problems encountered while running trace examples HOT 5
- Error when run xrtQook on Nvidia V100 GPU node
- How to add a caliber setting!
- Some doubts about the energy setting of X-ray sources HOT 18
- Absorbing wall along the beam propagation direction?
- Help needed! Clarification on the Meaning and Origin of eAxesFlat Values in 07_AnalyzerBent2D.py HOT 2
- Rocking curve theta and energy range converter HOT 3
- Problems generating code from materials HOT 7
- Calculating crystal reflectivity without constructing a beamline HOT 7
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
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.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from xrt.