Skip to content

Commit

Permalink
fix: translate quadrature
Browse files Browse the repository at this point in the history
  • Loading branch information
jokasimr committed Sep 17, 2024
1 parent df1bab0 commit 3a35667
Showing 1 changed file with 9 additions and 0 deletions.
9 changes: 9 additions & 0 deletions src/scippneutron/absorption/cylinder.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,10 @@ def beam_intersection(self, start_point, direction):
sc.scalar(0.0, unit=start_point.unit),
)

@property
def center(self):
return self.center_of_base + self.symmetry_line * self.height / 2

def _select_quadrature_points(self, kind):
if kind == 'expensive':
x, w = leggauss(max(round(10 * self.height / self.radius), 10))
Expand All @@ -56,9 +60,11 @@ def _select_quadrature_points(self, kind):

def quadrature(self, kind: Literal['expensive', 'medium', 'cheap']):
quad = self._select_quadrature_points(kind)
# Scale to size of cylinder
quad['x'] = quad['x'] * self.radius
quad['y'] = quad['y'] * self.radius
quad['z'] = quad['z'] * self.height / 2
quad['weights'] = quad['weights'] * (self.radius**2 * self.height / 2)
quad_points = sc.vectors(
dims=['quad'],
values=(
Expand All @@ -73,6 +79,9 @@ def quadrature(self, kind: Literal['expensive', 'medium', 'cheap']):
# We need to rotate the quadrature so the symmetry axis matches the cylinder.
u = sc.cross(sc.vector([0, 0, 1]), self.symmetry_line)
u *= sc.asin(sc.norm(u)) / sc.norm(u)
# By default the cylinder quadrature center is at the origin.
# We need to move it so the center matches the cylinder.
quad_points += self.center
return sc.spatial.rotations_from_rotvecs(u) * quad_points, sc.array(
dims=['quad'], values=quad['weights']
)
Expand Down

0 comments on commit 3a35667

Please sign in to comment.