Coordinate transformations

Coordinate transformations#

We often need to compute new coordinates from other coordinates; for example wavelength from time-of-flight or scattering angle from positions. Scipp has a flexible utility for this purpose: scipp.transform_coords.

Setup#

Consider a beamline with straight neutron beams (i.e., without guides, focusing optics, analyzers, etc.):

image

We want to compute the total length of the flight path Ltotal from positions of the detector, sample, and source. It is easy enough to write code that does this. However, consider the backscattering QENS beamline from the McStas session:

image

Here, we need to take the analyzer into account when computing Ltotal.

To illustrate, we begin with the example of straight beams and generate some test data. The exact contents of the data don’t matter here, but note that we store the various positions as coordinates.

import scipp as sc

# elastic_data
source_position = sc.vector([0.0, 0.0, -10.0], unit="m")
sample_position = sc.vector([0.0, 0.0, 0.0], unit="m")
position = sc.vectors(
    dims=["position"],
    values=[
        [0.0, 0.5, 1.0],
        [0.0, 1.0, 1.0],
        [0.0, 1.5, 1.0],
    ],
    unit="m",
)

elastic_data = sc.DataArray(
    sc.ones(sizes={"position": 3}),
    coords={
        "source_position": source_position,
        "sample_position": sample_position,
        "position": position,
    },
)
elastic_data
Show/Hide data repr Show/Hide attributes
scipp.DataArray (1.68 KB)
    • position: 3
    • position
      (position)
      vector3
      m
      [0. 0.5 1. ], [0. 1. 1.], [0. 1.5 1. ]
      Values:
      array([[0. , 0.5, 1. ],
       [0. , 1. , 1. ],
       [0. , 1.5, 1. ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -10.]
      Values:
      array([  0.,   0., -10.])
    • (position)
      float64
      𝟙
      1.0, 1.0, 1.0
      Values:
      array([1., 1., 1.])
# qens_data
source_position = sc.vector([0.0, 0.0, -10.0], unit="m")
sample_position = sc.vector([0.0, 0.0, 0.0], unit="m")
analyzer_position = sc.vector([0.0, 1.0, 1.0], unit="m")
position = sc.vectors(
    dims=["position"],
    values=[
        [0.0, 1.9, 0.0],
        [0.0, 2.0, 0.0],
        [0.0, 2.1, 0.0],
    ],
    unit="m",
)

qens_data = sc.DataArray(
    sc.ones(sizes={"position": 3}),
    coords={
        "source_position": source_position,
        "sample_position": sample_position,
        "analyzer_position": analyzer_position,
        "position": position,
    },
)
qens_data
Show/Hide data repr Show/Hide attributes
scipp.DataArray (2.00 KB)
    • position: 3
    • analyzer_position
      ()
      vector3
      m
      [0. 1. 1.]
      Values:
      array([0., 1., 1.])
    • position
      (position)
      vector3
      m
      [0. 1.9 0. ], [0. 2. 0.], [0. 2.1 0. ]
      Values:
      array([[0. , 1.9, 0. ],
       [0. , 2. , 0. ],
       [0. , 2.1, 0. ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -10.]
      Values:
      array([  0.,   0., -10.])
    • (position)
      float64
      𝟙
      1.0, 1.0, 1.0
      Values:
      array([1., 1., 1.])

The manual approach#

A straight-forward way of computing Ltotal from this is the following:

L1 = sc.norm(
    elastic_data.coords["sample_position"] - elastic_data.coords["source_position"]
)
L2 = sc.norm(elastic_data.coords["position"] - elastic_data.coords["sample_position"])
Ltotal = L1 + L2
Ltotal
Show/Hide data repr Show/Hide attributes
scipp.Variable (280 Bytes)
    • (position: 3)
      float64
      m
      11.118, 11.414, 11.803
      Values:
      array([11.11803399, 11.41421356, 11.80277564])

This uses vector arithmetic on the coordinates and scipp.norm to compute vector lengths.

However, if we now want to do this for the QENS experiment, we need to rewrite the entire procedure:

L1 = sc.norm(qens_data.coords["sample_position"] - qens_data.coords["source_position"])
L2 = sc.norm(
    qens_data.coords["position"] - qens_data.coords["analyzer_position"]
) + sc.norm(qens_data.coords["analyzer_position"] - qens_data.coords["sample_position"])
Ltotal = L1 + L2
Ltotal
Show/Hide data repr Show/Hide attributes
scipp.Variable (280 Bytes)
    • (position: 3)
      float64
      m
      12.760, 12.828, 12.901
      Values:
      array([12.75957597, 12.82842712, 12.90082044])

Using transform_coords#

Instead, we are going to use scipp.transform_coords.

First, we have to define functions to compute Ltotal and its components L1 and L2:

def straight_l1(source_position, sample_position):
    return sc.norm(sample_position - source_position)


def straight_l2(sample_position, position):
    return sc.norm(position - sample_position)


def l_total(L1, L2):
    return L1 + L2

We then store those functions in a dict. The dict-keys are names for the outputs of the functions.

graph = {"L1": straight_l1, "L2": straight_l2, "Ltotal": l_total}

This dict defines a graph that connects coordinates with functions that can compute them. We can visualize it with Scipp:

sc.show_graph(graph)
../_images/0ec226f2e427c1a646c467870c51985bc06d856c50dfd710677e5db9af9a2b6b.svg

Note how coordinates (white boxes) and functions (gray boxes) are connected. Scipp knows that, e.g., straight_l1 produces L1 because of the dict-key and uses source_position and sample_position as inputs because of the names of the function arguments.

We can now compute Ltotal by using the graph with transform_coords:

converted = elastic_data.transform_coords("Ltotal", graph=graph)
converted
Show/Hide data repr Show/Hide attributes
scipp.DataArray (2.48 KB)
    • Ltotal: 3
    • L1
      ()
      float64
      m
      10.0
      Values:
      array(10.)
    • L2
      (Ltotal)
      float64
      m
      1.118, 1.414, 1.803
      Values:
      array([1.11803399, 1.41421356, 1.80277564])
    • Ltotal
      (Ltotal)
      float64
      m
      11.118, 11.414, 11.803
      Values:
      array([11.11803399, 11.41421356, 11.80277564])
    • position
      (Ltotal)
      vector3
      m
      [0. 0.5 1. ], [0. 1. 1.], [0. 1.5 1. ]
      Values:
      array([[0. , 0.5, 1. ],
       [0. , 1. , 1. ],
       [0. , 1.5, 1. ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -10.]
      Values:
      array([  0.,   0., -10.])
    • (Ltotal)
      float64
      𝟙
      1.0, 1.0, 1.0
      Values:
      array([1., 1., 1.])

This did several things for us.

  • It computed Ltotal as we requested and stored it as a new coordinate.

  • It also computed L1 and L2 because those were needed for Ltotal.

  • It renamed the dimension from position to Ltotal because we consider the latter to have replaced the former.

Customizing the graph for QENS#

We can now adapt the above example to compute Ltotal for the QENS experiment. We need a new function that computes L2 the flight path length from sample to analyzer to detector.

image

def backscattering_l2(sample_position, analyzer_position, position):
    a = sc.norm(analyzer_position - sample_position)
    b = sc.norm(position - analyzer_position)
    return a + b

We can reuse the graph for the straight beamline and simply replace the function for L2:

graph["L2"] = backscattering_l2
sc.show_graph(graph)
../_images/56cf3e13c934b51618e3f59005416a68613c689a0fcfba1fe5bcc60608ef3155.svg 
converted = qens_data.transform_coords("Ltotal", graph=graph)
converted
Show/Hide data repr Show/Hide attributes
scipp.DataArray (2.81 KB)
    • Ltotal: 3
    • L1
      ()
      float64
      m
      10.0
      Values:
      array(10.)
    • L2
      (Ltotal)
      float64
      m
      2.760, 2.828, 2.901
      Values:
      array([2.75957597, 2.82842712, 2.90082044])
    • Ltotal
      (Ltotal)
      float64
      m
      12.760, 12.828, 12.901
      Values:
      array([12.75957597, 12.82842712, 12.90082044])
    • analyzer_position
      ()
      vector3
      m
      [0. 1. 1.]
      Values:
      array([0., 1., 1.])
    • position
      (Ltotal)
      vector3
      m
      [0. 1.9 0. ], [0. 2. 0.], [0. 2.1 0. ]
      Values:
      array([[0. , 1.9, 0. ],
       [0. , 2. , 0. ],
       [0. , 2.1, 0. ]])
    • sample_position
      ()
      vector3
      m
      [0. 0. 0.]
      Values:
      array([0., 0., 0.])
    • source_position
      ()
      vector3
      m
      [ 0. 0. -10.]
      Values:
      array([  0.,   0., -10.])
    • (Ltotal)
      float64
      𝟙
      1.0, 1.0, 1.0
      Values:
      array([1., 1., 1.])

The larger picture#

The examples shown above are fairly small and easy to see through. But in practice, coordinate transformations can involve more and more complicated steps. As an example, here is the default graph provided by ScippNeutron:

import scippneutron as scn

graph = scn.conversion.graph.beamline.beamline(scatter=True)
sc.show_graph(graph)
../_images/644bc63b98f2f4d0500363292daa5b957b42a508a8f127b350f71317f096c4aa.svg

It is similar to our own graph but involves additional intermediate results and can also be used to compute the scattering angle two_theta.

We can also add functions to compute the many more coordinates such as \(d\)-spacing, \(Q\), or hkl indices: (If you don’t know the syntax, simple read {**a, **b} as merging the two dicts a and b into a single dict.)

graph = {
    **scn.conversion.graph.beamline.beamline(scatter=True),
    **scn.conversion.graph.tof.elastic(start="tof"),
}
sc.show_graph(graph)
../_images/b27b0179037a95a452fc4b8c53656bc9981edf15c1fb8a50995e4d6befad5499.svg