[jl-wynen]

Fixes #730

This PR adds functions for

• basic handling of lattice parameters
• computing UB matrices
• goniometer definition and angle calculation

For now, this only supports spectroscopy because SXD uses a different definition of U. See in particular the docs of u_matrix_from_uv_vectors which describe this in terms of the orientation matrix $O$.

[jokasimr]

Is nan the only possible indication the matrix is singular, or might we see inf as well? In that case np.all(np.isfinite(...)) might be the better check.

[jl-wynen]

Good question. The behaviour is undefined: https://libeigen.gitlab.io/eigen/docs-nightly/classEigen_1_1MatrixBase.html#a7712eb69e8ea3c8f7b8da1c44dbdeebf So it could return NaN, Inf, 4, or anything else.

We only use this in one place now (inverting $G$) and it is not even clear if we need that. I could just remove it for now?

[jokasimr]

Do what you think is best 👍

[jokasimr]

I'm not sure I understood this. Are these the conditions determining $U$?

1. $U$ is a rotation:

$$ U \in SO(3) $$

2. $U$ maps the provided index to the corresponding momentum transfer:

$$ q_s = 2\pi\ U\ B\ h $$

3. $e_x$ in the sample system is in the scattering plane:

$$ \begin{bmatrix}1 & 0 & 0\end{bmatrix} U s_n = 0 $$

4. $e_y$ in the sample system is "aligned" with the scattering plane normal:

$$ \begin{bmatrix}0 & 1 & 0\end{bmatrix} U s_n > 0 $$

[jl-wynen]

I think I improved the wording. Can you check if it is clearer now?

[jokasimr]

In the formula Q_nu^T = V^T * (UB)^T it looks like the solution should be solve(V.T, q_nu.T) but instead here it says solve(V, q_nu).
Are q_nu.values and peaks.hkl.values already transposed?

[jl-wynen]

Yes. In matrix notation, $V$ and $UB$ contain three column vectors. But we have variables with dtype=vector3 and their .values has the vectors as rows, i.e., transposed compared to the matrix notation.

[jokasimr]

align_vectors allows passing a weight argument where one weight can be infinite. That seems to solve a similar (the same?) problem as here, where a rotation is under determined but we use a minimization condition to make it determined.

Would it be possible to use that mechanism here instead of implementing the constrained minimization manually?

[jl-wynen]

Can you elaborate? How would this resolve the ambiguity?

[jokasimr]

If I understand correctly this function determines a rotation $R$ such that the peak at $hkl$ is mapped exactly to the expected $q$ vector in the lab system, at the given scattering angle. But that rotation is ambiguous, so to resolve the ambiguity we select the rotation that minimizes the "tilt" of the sample. That is:

$$ R = arg \ max_R \quad \vec{e}_y \cdot R \vec{n} $$ $$ s.t. \quad \vec{q}(\theta, e_i, e_f) = RUB\ \vec{h} $$

where $\vec{n}$ is the orientation of the sample.

I think the above optimization problem can be solved directly by scipy.Rotation.align_vectors using something like:

R = align_vectors(
    [q(theta, ei, ef), e_y]
    [ub * hkl, n],
    weights=[np.inf, 1.0]  # Note the infinite weight.
)

So there might be opportunity for a simplification of this function.

[jl-wynen]

This doesn't work. With the example in the docs, we get these angles: (left is old, right is your idea)

+x : +0.0000 deg   |  -90.0009 deg
-z : -5.8930 deg   |  -0.0085 deg
+y : -46.9236 deg  |  -41.0220 deg

And rotating the vector [0, 1, 0] by these rotations gives

old: [-0.07012154 0.99471537 -0.07499537]
new: [6.56348769e-01 -1.53888876e-05 -7.54457615e-01]

So the new method does not try to minimise the tilt.

I tried with a weight smaller than Inf but then the result does not converge. I guess because it needs 3 vectors to converge reliably.

[jokasimr]

But it shouldn't need three vectors to converge reliably? It should only need two 🤔 I mean, each vector contributes two equations and the rotation has three degrees of freedom.

So I guess the minimization problem that I wrote down is not the one solved by this function?

[jl-wynen]

I accidentally used the wrong vector for e_y. When I use the correct one, it works. I tested with a bunch of different hkl vectors and the result is always within 1e-16 of the current implementation. I'll update the code. Good find!

[jl-wynen]

On second thought, this doesn't work when hkl is parallel to the plane normal. The existing impl has a case to handle that, but not the new one. And they produce different results:

hkl = sc.vector([1, 1, 0])

# Make a plane with normal || ub * hkl
n = ub1 * hkl
xx = 0.12
zz = -n.fields.x * xx / n.fields.z
v1 = sc.vector([xx, 0, zz.value])
plane = ScatteringPlaneSampleHolder.from_normal_and_reference(normal=n, reference=v1)

old = alignment.target_rotation_for_reflection(
    hkl=hkl,
    two_theta=known_two_theta2,
    ei=ei,
    ef=ei,  # We only care about elastic scattering for locating peaks.
    ub=ub1,
    scattering_plane=plane,
)
new = new_impl(
    hkl=hkl,
    two_theta=known_two_theta2,
    ei=ei,
    ef=ei,  # We only care about elastic scattering for locating peaks.
    ub=ub1,
    scattering_plane=plane,
)
print((old * sc.vector([0, 1, 0])).value)
print((new * sc.vector([0, 1, 0])).value)
print((old * sc.vector([0, 0, 1])).value)
print((new * sc.vector([0, 0, 1])).value)

[ 0.05953926 -0.9961941   0.06365843]
[ 0.08710432 -0.99547886  0.03787704]
[ 7.30263991e-01 -9.88036077e-06 -6.83165063e-01]
[ 0.7297399   0.03787704 -0.68267489]

I'd rather stick with what we have because that was based on TAVI which has been used in practice and probably produces better results in these edge cases.

[jokasimr]

Let's stick with what we have in that case!

[jokasimr]

Suggested change

	q_hkl2p = b * scattering_plane.vector1 * np.sign(c2.value)
	q_hkl2p = -b * scattering_plane.vector1 * np.sign(c2.value)

because otherwise it seems, if v1=[0, 0, 1] and v2=[1, 0, 0] and q = [1, 0, 1], then U will map [0, 1, 0] in the crystal frame to [0, -1, 0] in the sample frame.

[jl-wynen]

I don't understand what you mean. We should talk about this in person.