Incorrect calculation of exit angle for Q vector

[P-Mousley]

I think that the calculation performed when requesting Q units for 2d map calculations appears to be only valid when x=0, which does not work if you want to move the detector using rot2. The function eq_exitangle that has the issue is here

the current version is simply

numpy.arctan2(y, z) - incident_angle

which does not take into account the effect of x, and also assumes that the exit vector is aligned with the incident angle

a possible correction is as follows:

def eq_exitangle(x, y, z, wavelength=None, incident_angle=0.0, tilt_angle=0.0, sample_orientation=1):
    """Calculates the vertical exit scattering angle (relative to incident angle), used for GI/Fiber diffraction
    :param x: horizontal position, towards the center of the ring, from sample position
    :param y: vertical position, to the roof, from sample position
    :param z: distance from sample along the beam
    :param wavelength: in meter
    :return: exit angle for the position xyz in radians
    """
    #add in correction to do proper angle calculation - first need to get vector length]
    #then find height in y to use sin rule for angle to the y=0 plane
    fulllength=numpy.linalg.norm((x,y,z))
    #also need proper calculation of adjustment to height dependent on tilt_angle
    th1=numpy.arctan(z/y)
    height=y*(numpy.cos(th1+tilt_angle)/numpy.cos(th1))
    #use new values to calulate the radian angle to the xz plane at y=0
    rad_to_xz=numpy.arcsin(height/fulllength)
    return rad_to_xz

Let me know if i have understood the function correctly, and if the proposed correction is suitable.

[kif]

Hi Philip, Thanks for the report, I will let Edgar deal with this since I am not the most competent on this topic. Cheers, Jerome

[EdgarGF93]

Hi, Philip. Thanks a lot for the issue. Indeed, there is an inconsistency with the exit_angle equation. However, to be consistent with the way I build the q equations, the correct equation is:

 numpy.arctan2(y, numpy.sqrt(z ** 2 + x ** 2))

I use this exit_angle before applying the incident_angle and tilt_angle rotation matrix, so this exit_angle is a magnitude relative to the direct beam axis: it's the angle between kf and kf_xy

So, you're right that it depends on x, but it's a magnitude of the pixel position (only dependent of the poni and independent of the sample rotations). Again, this is my 'exit_angle' that I define like this, and then I apply the two rotation matrix to transform the q vector into the sample reference.

These days, I will include these corrections, also I will clarify the API with the equations, which are quite bulky.

[EdgarGF93]

#2297

[EdgarGF93]

Although I understand that the grazing incidence community calls exit_angle a magnitude relative to the horizon of the thin film sample, so it's a magnitude of the sample reference frame that, of course, depends on both incident_angle and tilt_angle. For my approach, what I call now vertical and horizontal exit angles are actually the components of the classic scattering_angle (2theta) from the powder diffraction point of view, which is I think the more "pyFAInic" way to implement this. It would be happy anyway to include this grazing_incidence exit_angle unit in the API.

[P-Mousley]

Thank you for the explanation - I think the new vector component calculations make sense to me. If i have understood correctly the rot1,rot2,rot3 from the PONI file are used to update the XYZ in the labframe for the exit vector using equations defined in pyFAI/src/pyFAI/ext/_geometry.pyx . The updated XYZ is then passed to your updated equations in units.py to calculated the Q vector components, taking account of incident_angle and tilt_angle

i still have two questions:

1. Am i correct in assuming this mean that to provide the details for incident_angle and tilt_angle you can no longer use

unit_qip = "qip_nm^-1"
unit_qoop = "qoop_nm^-1"

instead you need to use the grazing incidence fiber units e.g.

unit_qip = get_unit_fiber(name="qip_nm^-1", incident_angle=0.2, tilt_angle=0.4, sample_orientation=1)
unit_qoop = get_unit_fiber(name="qoop_nm^-1", incident_angle=0.2, tilt_angle=0.4, sample_orientation=1)

2. Is there a specific reason why rotation 1 (around the axis towards the celing) needs to be done before rotation 2 (around the axis pointing towards the ring)? Because this makes it difficult to simulate a diffractometer where the out of plane angle (delta) is stacked ontop of the inplane angle (gamma). I think switching the order of how the angles are applied within the conversion equations would make it easier to translate from a (gamma, dellta) value, to a (rot1,rot2) value. However i am not sure if this would have larger consequences for the pyFAI codebase.

[kif]

Hi Philip, I am not a fiber guy so I can only answer for your second point. I wanted the last rotation to be free along Debye-Scherrer rings. There was no specific reason for the order of the rotation around vertical and horizontal axis of the detector and honestly, when I started pyFAI, I did not consider it would ever be used in conjunction with a diffractometer since everybody was happy with a detector mounted normal to the beam. Now pyFAI is used by thousands of users and PONI-files are distributed to all our users. Changing something so central is just not an option. But I can try to help in the math/trigonometry for converting geometries ... and this is something we would like to address for pyFAI2 (pluggable geometries)

[EdgarGF93]

#2308

[P-Mousley]

Hello, I have had a look at the conversion equations to get a rot1,rot2 value pair for a set value of gamma and delta. However it seems impossible to get a direct conversion because pyFAI applies the same rot1 and rot2 values to the whole image. When comparing the two methods of positioning the detector (either rot1 then rot2, or delta then gamma) there is always a mismatch with the image vectors even if the central pixel overlaps perfectly. the attached plot shows this mismatch:

• the blue square the detector at all angels at 0
• the orange square for rotation of delta then gamma
• the green square for rotation of rot1 then rot2.

Screencast from 29-10-24 15:41:43.webm

Do you have any more details of what other users have been doing if they want to use pyFAI for images taken with diffractometers where the delta motion is stacked ontop of the gamma motion?

[kif]

I believe that if the two rotation are actually inverted (rot1+rot2 vs gamma+delta), the solution is to use an Euler transformation library like the one from Christoph Gohlke which is able to match any of many set of Euler angle. https://pypi.org/project/transformations/ There is a snapshot of this library in pyFAI.third_party.transformations which should be usable for this purpose.

[mjdiff]

Hello @P-Mousley , could you share with us how you calculate coordinates in your script?

Hello, I have had a look at the conversion equations to get a rot1,rot2 value pair for a set value of gamma and delta. However it seems impossible to get a direct conversion because pyFAI applies the same rot1 and rot2 values to the whole image. When comparing the two methods of positioning the detector (either rot1 then rot2, or delta then gamma) there is always a mismatch with the image vectors even if the central pixel overlaps perfectly. the attached plot shows this mismatch:

* the blue square the detector at all angels at 0

* the orange square for rotation of delta then gamma

* the green square for rotation of rot1 then rot2.

Screencast.from.29-10-24.15.41.43.webm

Do you have any more details of what other users have been doing if they want to use pyFAI for images taken with diffractometers where the delta motion is stacked ontop of the gamma motion?

[P-Mousley]

Hello @P-Mousley , could you share with us how you calculate coordinates in your script?

Hello, I have had a look at the conversion equations to get a rot1,rot2 value pair for a set value of gamma and delta. However it seems impossible to get a direct conversion because pyFAI applies the same rot1 and rot2 values to the whole image. When comparing the two methods of positioning the detector (either rot1 then rot2, or delta then gamma) there is always a mismatch with the image vectors even if the central pixel overlaps perfectly. the attached plot shows this mismatch:

* the blue square the detector at all angels at 0

* the orange square for rotation of delta then gamma

* the green square for rotation of rot1 then rot2.

Screencast.from.29-10-24.15.41.43.webm Do you have any more details of what other users have been doing if they want to use pyFAI for images taken with diffractometers where the delta motion is stacked ontop of the gamma motion?

Hello @mjdiff ,

To calculate the co-ordinates in the script I used simple rotation matrices from scipy around the x,y,or z axis.

Expand code used for rotation calculation

import matplotlib.pyplot as plt
from scipy.spatial.transform import Rotation as R
import transformations as tf
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as wg
%matplotlib qt
ax = plt.axes([0,0,0.95,0.95],projection ='3d')


def drawline(vec1,vec2,axs,col):
    axs.plot([vec1[0],vec2[0]],[vec1[1],vec2[1]],[vec1[2],vec2[2]],color = col)

def drawdet(veclist,ax,col):
    drawline(veclist[0],veclist[1],ax,col)
    drawline(veclist[1],veclist[2],ax,col)
    drawline(veclist[2],veclist[3],ax,col)
    drawline(veclist[3],veclist[0],ax,col)
    drawline([0,0,0],veclist[4],ax,col)

def drawplots(gamma,delta,rot1,rot2,rot3):
    ax.cla()
    # Set the axis labels
    ax.set_xlabel('z')
    ax.set_ylabel('x')
    ax.set_zlabel('y')
    rotdelta=R.from_euler('y', -delta, degrees=True)
    rotgamma=R.from_euler('z',gamma,degrees=True)
    rot1mat=R.from_euler('y', -rot2, degrees=True)
    rot2mat=R.from_euler('z',rot1,degrees=True)
    rot3mat=R.from_euler('x',rot3,degrees=True)

    veclist=[[4,0.5,0.5],[4,0.5,-0.5],[4,-0.5,-0.5],[4,-0.5,0.5],[4,0,0]]
    totalrot=rotgamma*rotdelta
    totalrot2=rot3mat*rot2mat*rot1mat
    
    rotvecs1=[]
    rotvecs2=[]
    for vec in veclist:
        rotvecs1.append(totalrot.apply(vec))
        rotvecs2.append(totalrot2.apply(vec))


    drawdet(veclist,ax,'blue')
    drawdet(rotvecs1,ax,'orange')
    drawdet(rotvecs2,ax,'green')
    
    ax.set_zlim(-4,4)
    ax.set_ylim(-4,4)
    ax.set_xlim(-4,4)
gammaslide=wg.FloatSlider(min=-10, max=90, step=0.001, value=0)
deltaslide=wg.FloatSlider(min=-10, max=90, step=0.001, value=0)
rot1slide=wg.FloatSlider(min=-90, max=90, step=0.001, value=0)
rot2slide=wg.FloatSlider(min=-90, max=90, step=0.001, value=0)
rot3slide=wg.FloatSlider(min=-90, max=90, step=0.001, value=0)
interact(drawplots,axs=ax,gamma=gammaslide,delta=deltaslide,rot1=rot1slide,rot2=rot2slide,rot3=rot3slide)

The transformations library @kif recommend earlier has allowed me to calculate the corresponding values for rot1,rot2,rot3 which are equivalent to a delta +gamma pair of diffractometer angles - see code example below.

Expand code used for converting delta,gamma into rot1,rot2,rot3

    def gamdel2rots(self,gamma,delta):
        """
        

        Parameters
        ----------
        gamma : float
            angle rotation of gamma diffractometer circle in degrees.
        delta : float
            angle rotation of delta diffractometer circle in degrees.

        Returns
        -------
        rots : list of rotations rot1,rot2,rot3 in radians to be using by pyFAI.

        """
        rotdelta=R.from_euler('y', -delta, degrees=True)
        rotgamma=R.from_euler('z',gamma,degrees=True)
        totalrot=rotgamma*rotdelta
        fullrot=np.identity(4)
        fullrot[0:3,0:3]=totalrot.as_matrix()
        vals=tf.euler_from_matrix(fullrot,'rxyz')
        rots=vals[2],-vals[1],vals[0]
        return rots

[kif]

Great ... Thanks for your work. Can we reuse it in the documentation ?

[P-Mousley]

Great ... Thanks for your work. Can we reuse it in the documentation ?

Yes, happy for you to use my code

[Anas321]

Hello All. So, is this issue closed now? Was it corrected?

[kif]

I believe this is mostly integrated in the latest version of pyFAI but maybe some documentation remains to be updated, thus the unclosed issue. Edgar could comment on it but he is currently off.

[Anas321]

@kif Ok, thanks so much for your response.

I also have another, different, question. If I want to look at the equations used to construct the q vector, is this the right place to look in: https://github.com/silx-kit/pyFAI/blob/main/src/pyFAI/units.py

I am trying to compare what I am getting from pyFAI with my manual calculations using equation (1) in Hexemer's paper: https://journals.iucr.org/m/issues/2015/01/00/ed5003/ed5003.pdf

However, I am finding some discrepancies in the vertical direction of the q vector, and I just want to confirm that I am looking at the right place in pyFAI's source code. I am still investigating this. There might be a mistake on my side.

Appreciate your help. Thanks!

[kif]

Units.py is the right place to look for. I will also read the paper and check ... but I am not a specialist of surface diffraction.

[Anas321]

@kif Thanks so much for your help and prompt response.