Skip to content

Counter rotating choppers #44

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
nvaytet wants to merge 5 commits into from
Closed

Counter rotating choppers #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1b17177
add support for counter-rotating choppers via negative frequency
nvaytet Jun 9, 2023
7998929
Merge branch 'main' into counter-rotating-choppers
nvaytet Jun 9, 2023
588d571
fix unit tests
nvaytet Aug 7, 2023
4260d20
update docs with negative frequencies
nvaytet Aug 7, 2023
996c476
flake8
nvaytet Aug 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 9 additions & 3 deletions docs/components.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,13 @@
"\n",
"Next, we add a chopper in the beamline,\n",
"with a frequency, phase, distance from source,\n",
"and a set of open and close angles for the cutouts in the rotating disk."
"and a set of open and close angles for the cutouts in the rotating disk.\n",
"\n",
"Note that a chopper rotating anti-clockwise is given a positive frequency,\n",
"while a chopper rotating clockwise is given a negative frequency.\n",
Comment on lines +140 to +141
Copy link
Copy Markdown
Member

@SimonHeybrock SimonHeybrock Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this convention used?

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the convention from Nexus, but I don't like it either.
McStas also uses the sign of the frequency, and theirs is the opposite (positive for clockwise).

But regarding your comment below, yes I think I prefer using positive values for frequency and then a direction parameter.

"\n",
"If you define your window cutouts in the **anti-clockwise** direction,\n",
"your chopper should typically rotate in the **clockwise** direction so that the window with the lowest cutout angles is the first one to pass in front of the beam."
]
},
{
Expand All @@ -146,7 +152,7 @@
"outputs": [],
"source": [
"chopper1 = tof.Chopper(\n",
" frequency=10.0 * Hz,\n",
" frequency=-10.0 * Hz, # Negative frequency for clockwise rotation\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[30.0, 50.0],\n",
Expand Down Expand Up @@ -246,7 +252,7 @@
"outputs": [],
"source": [
"chopper2 = tof.Chopper(\n",
" frequency=5.0 * Hz,\n",
" frequency=-5.0 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[30.0, 40.0, 55.0, 70.0],\n",
Expand Down
12 changes: 6 additions & 6 deletions docs/multiple-pulses.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@
"source": [
"choppers = [\n",
" tof.Chopper(\n",
" frequency=70.0 * Hz,\n",
" frequency=-70.0 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[98.71, 155.49, 208.26, 257.32, 302.91, 345.3],\n",
Expand All @@ -91,7 +91,7 @@
" name=\"WFM1\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=70 * Hz,\n",
" frequency=-70 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[80.04, 141.1, 197.88, 250.67, 299.73, 345.0],\n",
Expand All @@ -107,7 +107,7 @@
" name=\"WFM2\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=56 * Hz,\n",
" frequency=-56 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[74.6, 139.6, 194.3, 245.3, 294.8, 347.2],\n",
Expand All @@ -123,7 +123,7 @@
" name=\"Frame-overlap 1\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=28 * Hz,\n",
" frequency=-28 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[98.0, 154.0, 206.8, 254.0, 299.0, 344.65],\n",
Expand Down Expand Up @@ -175,7 +175,7 @@
"metadata": {},
"outputs": [],
"source": [
"res.plot(max_rays=5000)"
"res.plot(max_rays=10000)"
]
},
{
Expand Down Expand Up @@ -216,7 +216,7 @@
"outputs": [],
"source": [
"pol = tof.Chopper(\n",
" frequency=14 * Hz,\n",
" frequency=-14 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[50.0],\n",
Expand Down
10 changes: 5 additions & 5 deletions docs/short-example.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@
"source": [
"choppers = [\n",
" tof.Chopper(\n",
" frequency=70.0 * Hz,\n",
" frequency=-70.0 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[98.71, 155.49, 208.26, 257.32, 302.91, 345.3],\n",
Expand All @@ -93,7 +93,7 @@
" name=\"WFM1\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=70 * Hz,\n",
" frequency=-70 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[80.04, 141.1, 197.88, 250.67, 299.73, 345.0],\n",
Expand All @@ -109,7 +109,7 @@
" name=\"WFM2\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=56 * Hz,\n",
" frequency=-56 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[74.6, 139.6, 194.3, 245.3, 294.8, 347.2],\n",
Expand All @@ -125,7 +125,7 @@
" name=\"Frame-overlap 1\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=28 * Hz,\n",
" frequency=-28 * Hz,\n",
" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[98.0, 154.0, 206.8, 254.0, 299.0, 344.65],\n",
Expand All @@ -141,7 +141,7 @@
" name=\"Frame-overlap 2\",\n",
" ),\n",
" tof.Chopper(\n",
" frequency=7 * Hz,\n",
" frequency=-7 * Hz,\n",
Comment on lines 143 to +144
Copy link
Copy Markdown
Member

@SimonHeybrock SimonHeybrock Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if relying on a sign convention will lead to too many bugs? Have you considered splitting this into frequency (positive) and direction: Literal['clockwise', 'anticlockwise']?

" open=sc.array(\n",
" dims=['cutout'],\n",
" values=[30.0],\n",
Expand Down
45 changes: 37 additions & 8 deletions src/tof/chopper.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
from dataclasses import dataclass
from typing import Optional, Tuple

import numpy as np
import scipp as sc

from .reading import ComponentReading, ReadingField
Expand All @@ -18,15 +19,19 @@ class Chopper:
Parameters
----------
frequency:
The frequency of the chopper.
The frequency of the chopper. Positive values indicate anti-clockwise rotation,
negative values indicate clockwise rotation.
open:
The opening angles of the chopper cutouts.
close:
The closing angles of the chopper cutouts.
distance:
The distance from the source to the chopper.
phase:
The phase of the chopper.
The phase of the chopper. Note that the phase is applied opposite to the
direction of rotation. For example, if the chopper rotates clockwise, a
phase of 10 degrees will shift all the openings by 10 degrees in the
anti-clockwise direction.
name:
The name of the chopper.
"""
Expand Down Expand Up @@ -76,19 +81,43 @@ def open_close_times(
"""
if unit is None:
unit = time_limit.unit
nrot = max(int(sc.ceil((time_limit * self.frequency).to(unit='')).value), 1)
nrot = max(
int(sc.ceil((time_limit * abs(self.frequency)).to(unit='')).value), 1
)
# Start at -1 to catch early openings in case the phase or opening angles are
# large
phases = sc.arange(uuid.uuid4().hex, -1, nrot) * two_pi + self.phase.to(
unit='rad'
)
) * (np.sign(self.frequency.value) * -1.0)
# Note that the order is important here: we need (phases + open/close) to get
# the correct dimension order when we flatten below.
open_times = (phases + self.open.to(unit='rad', copy=False)) / self.omega
close_times = (phases + self.close.to(unit='rad', copy=False)) / self.omega
open_times = (phases + self.open.to(unit='rad', copy=False)).flatten(
to=self.open.dim
)
close_times = (phases + self.close.to(unit='rad', copy=False)).flatten(
to=self.close.dim
)
# If the chopper is rotating anti-clockwise, we mirror the openings because the
# first cutout will be the last to open.
if self.frequency.value > 0:
open_times, close_times = (
sc.array(
dims=close_times.dims,
values=(two_pi - close_times).values[::-1],
unit=close_times.unit,
),
sc.array(
dims=open_times.dims,
values=(two_pi - open_times).values[::-1],
unit=open_times.unit,
),
)
Comment on lines +100 to +114
Copy link
Copy Markdown
Member

@SimonHeybrock SimonHeybrock Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The interface and need to handle this here feels very error-prone. Would it be possible to always specify openings in the order they appear, with positive angles? Any code for converting from different convention could be moved into helpers, such as Angles.from_nexus(**nexus_chopper_angles).

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice suggestion, but wouldn't this mean that the direction parameter would then be obsolete? Because we would always have the first opening crossing the beam first, regardless of the direction of rotation?

So I am wondering if it could me more error-prone in some cases?

Copy link
Copy Markdown
Member

@SimonHeybrock SimonHeybrock Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure how different conventions for different things (rotation, phase, angles) interact, so this was just an idea. I don't mind either solution, but I thought that direction alone might not be sufficient, because there may be other conventions as well?

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Aug 7, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, now that I think a bit more about this, I am wondering if we should do what you suggest, and not have a sign on the frequency, nor a direction parameter...

In the case of the dream choppers, there are two counter-rotating choppers, where one is above the beam and the other one is below. There is even a third chopper which is on the side:
Screenshot at 2023-08-07 17-51-31

  • NeXus defines a TDC on the chopper, and an offset between the TDC and the position of the beam.
  • McStas deals with this by defining the centre of the beam as the zero angle, meaning that the chopper above the beam is rotated 180deg compared to the one below.

I think defining things in terms of "the first window that is listed is the first one that will pass in front of the beam" makes sense. If a user has angles from a chopper diagram, and then they want to place that chopper on the side of the beam, they just have to either add 90deg to all their angles, or add 90deg to their phase.

To define the opening angles for a counter-rotating chopper, you have to reverse the order from what is in the NeXus.
I like the idea of having converters in separate functions.

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Aug 7, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this mean we have to change nothing to the code but instead update the documentation to explain how to make a counter-rotating chopper?

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Aug 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this mean we have to change nothing to the code but instead update the documentation to explain how to make a counter-rotating chopper?

I think we probably still need a direction parameter. Users typically have a diagram giving the opening angles:
Screenshot at 2023-08-07 18-46-27

It will not always be from_nexus or from_mcstas, sometimes you just want to check the values given in a paper. I think it will be useful to be able to set direction='anticlockwise'.

Copy link
Copy Markdown
Member Author

@nvaytet nvaytet Aug 7, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

About Angles.from_nexus(**nexus_chopper_angles), wouldn't we try and load everything from NeXus, not just the angles but the frequency, the phase etc?

For McStas, we cannot parse the instrument file, so I am not sure what to do there.

I am asking because I thought about either having Angles.counter_rotating(**angles), or a method on the chopper to flip the rotation, e.g.

chopper = Chopper(...)
chopper.counter_rotate()  # would modify in-place the open and close angles

The problem with the latter is that you don't necessarily remember how many times you called counter_rotate() and if you are doing stuff in the notebook, you may accidentally flip your chopper again by executing a cell once too many, so it would probably be better to have is set when constructing the chopper.

How about using classmethods?

  • Chopper.anti_clockwise(...)
  • Chopper.from_nexus(...)

Copy link
Copy Markdown
Member

@SimonHeybrock SimonHeybrock Aug 8, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 on the classmethods, this is more or less what I had in mind.

Yes, the idea was to not change this code here, but instead separate those concerns into independent helpers.

abs_omg = abs(self.omega)
open_times /= abs_omg
close_times /= abs_omg
return (
open_times.flatten(to=self.open.dim).to(unit=unit, copy=False),
close_times.flatten(to=self.close.dim).to(unit=unit, copy=False),
open_times.to(unit=unit, copy=False),
close_times.to(unit=unit, copy=False),
)

def __repr__(self) -> str:
Expand Down
91 changes: 82 additions & 9 deletions tests/chopper_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,22 @@ def test_angular_speed():
assert chopper.omega == two_pi * f


def test_angular_speed_negative():
f = -8.0 * Hz
chopper = tof.Chopper(
frequency=f,
open=0.0 * deg,
close=10.0 * deg,
phase=0.0 * deg,
distance=5.0 * meter,
)
assert chopper.omega == two_pi * f


def test_open_close_times_one_rotation():
f = 10.0 * Hz
chopper = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=sc.array(dims=['cutout'], values=[10.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[20.0], unit='deg'),
phase=0.0 * deg,
Expand All @@ -49,7 +61,7 @@ def test_open_close_times_one_rotation():
def test_open_close_times_three_rotations():
f = 10.0 * Hz
chopper = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=sc.array(dims=['cutout'], values=[10.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[20.0], unit='deg'),
phase=0.0 * deg,
Expand Down Expand Up @@ -79,7 +91,7 @@ def test_open_close_times_three_rotations():
def test_open_close_angles_scalars_converted_to_arrays():
f = 10.0 * Hz
chopper = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=10.0 * deg,
close=20.0 * deg,
phase=0.0 * deg,
Expand All @@ -93,23 +105,27 @@ def test_open_close_angles_scalars_converted_to_arrays():
def test_phase():
f = 10.0 * Hz
chopper1 = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=sc.array(dims=['cutout'], values=[10.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[20.0], unit='deg'),
phase=0.0 * deg,
distance=10.0 * meter,
)
topen1, tclose1 = chopper1.open_close_times(0.0 * sec)
chopper2 = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=sc.array(dims=['cutout'], values=[10.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[20.0], unit='deg'),
phase=30.0 * deg,
distance=10.0 * meter,
)
topen2, tclose2 = chopper2.open_close_times(0.0 * sec)
assert sc.allclose(topen2, topen1 + (30.0 * deg).to(unit='rad') / chopper2.omega)
assert sc.allclose(tclose2, tclose1 + (30.0 * deg).to(unit='rad') / chopper2.omega)
assert sc.allclose(
topen2, topen1 + (30.0 * deg).to(unit='rad') / abs(chopper2.omega)
)
assert sc.allclose(
tclose2, tclose1 + (30.0 * deg).to(unit='rad') / abs(chopper2.omega)
)


def test_phase_int():
Expand All @@ -118,14 +134,14 @@ def test_phase_int():
cl = sc.array(dims=['cutout'], values=[20.0], unit='deg')
d = 10.0 * meter
chopper1 = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=op,
close=cl,
phase=30.0 * deg,
distance=d,
)
chopper2 = tof.Chopper(
frequency=f,
frequency=-f, # negative frequency for clockwise rotation
open=op,
close=cl,
phase=30 * deg,
Expand All @@ -135,3 +151,60 @@ def test_phase_int():
topen2, tclose2 = chopper2.open_close_times(0.0 * sec)
assert sc.allclose(topen1, topen2)
assert sc.allclose(tclose1, tclose2)


def test_open_close_times_counter_rotation():
f = 10.0 * Hz
d = 10.0 * meter
ph = 0.0 * deg
chopper1 = tof.Chopper(
frequency=-f,
open=sc.array(dims=['cutout'], values=[10.0, 90.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[20.0, 130.0], unit='deg'),
phase=ph,
distance=d,
)
chopper2 = tof.Chopper(
frequency=f,
open=sc.array(
dims=['cutout'], values=[360.0 - 130.0, 360.0 - 20.0], unit='deg'
),
close=sc.array(
dims=['cutout'], values=[360.0 - 90.0, 360.0 - 10.0], unit='deg'
),
phase=ph,
distance=d,
)

topen1, tclose1 = chopper1.open_close_times(0.2 * sec)
topen2, tclose2 = chopper2.open_close_times(0.0 * sec)
# Note that the first chopper will have one more rotation before t=0, so we slice
# out the first two open/close times
assert sc.allclose(topen1[2:], topen2)
assert sc.allclose(tclose1[2:], tclose2)


def test_open_close_times_counter_rotation_with_phase():
f = 10.0 * Hz
chopper1 = tof.Chopper(
frequency=f,
open=sc.array(dims=['cutout'], values=[80.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[90.0], unit='deg'),
phase=0.0 * deg,
distance=10.0 * meter,
)
topen1, tclose1 = chopper1.open_close_times(0.0 * sec)
chopper2 = tof.Chopper(
frequency=f,
open=sc.array(dims=['cutout'], values=[80.0], unit='deg'),
close=sc.array(dims=['cutout'], values=[90.0], unit='deg'),
phase=30.0 * deg,
distance=10.0 * meter,
)
topen2, tclose2 = chopper2.open_close_times(0.0 * sec)
assert sc.allclose(
topen2, topen1 + (30.0 * deg).to(unit='rad') / abs(chopper2.omega)
)
assert sc.allclose(
tclose2, tclose1 + (30.0 * deg).to(unit='rad') / abs(chopper2.omega)
)
4 changes: 2 additions & 2 deletions tests/common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@


def make_chopper(topen, tclose, f, phase, distance, name):
aopen = two_pi * sc.concat(topen, dim='cutout').to(unit='s') * f
aclose = two_pi * sc.concat(tclose, dim='cutout').to(unit='s') * f
aopen = two_pi * sc.concat(topen, dim='cutout').to(unit='s') * (-f)
aclose = two_pi * sc.concat(tclose, dim='cutout').to(unit='s') * (-f)
return tof.Chopper(
frequency=f,
open=aopen,
Expand Down
Loading