clatlan / cdiutils

A python package to help Coherent Diffraction Imaging (CDI) practitioners in their analysis.
MIT License
12 stars 7 forks source link

[FEATURE] tilt computation for a Bragg not along one of the axes #27

Open ewbellec opened 9 months ago

ewbellec commented 9 months ago

In cases where the Bragg is not along one of the array axes, the tilt is not directly given by the gradient components along the other axes. In that case, you need to find 2 vectors perpendicular to the Bragg (here I use orthogonal_vectors function on the Bragg wavevector qcen) and project the gradient (grad) on these vectors in compute_tilt function. These are just to give some ideas. Use whatever you want.

def orthogonal_vectors(qcen,
                      check=False):
    '''
    :qcen: Bragg vector (center of mass of the BCDI peak)
    '''
    e_bragg = qcen / np.linalg.norm(qcen)

    guess = np.zeros(3)
    guess[np.argmin(qcen)] = 1

    e1 = np.cross(guess, e_bragg)
    e1 = e1 / np.linalg.norm(e1)

    e2 = np.cross(e_bragg, e1)
    e2 = e2 / np.linalg.norm(e2)

    if check:
        print('e_bragg : ', e_bragg)
        print('e1 : ', e1)
        print('e2 : ', e2)

        print('all values below should be 1')
        print('check 1', np.dot(e_bragg, e_bragg))
        print('check 2', np.dot(e1, e1))
        print('check 3', np.dot(e2, e2))

        print('\nall values below should be 0')
        print('check 4', np.dot(e_bragg, e1))
        print('check 5', np.dot(e_bragg, e2))
        print('check 6', np.dot(e1, e2))
    return e_bragg, e1, e2

def compute_tilt(grad, qcen,
                 polar_representation=True,
                 check=False):

    # Get 2 vectors perpendicular to the Bragg
    e_bragg, e1, e2 = orthogonal_vectors(qcen, check=check)

    # components of  the tilt along e1 and e2
    tilt_comp1 =  np.sum([e1[n] * grad[n] for n in range(len(qcen))], axis=0)
    tilt_comp2 =  np.sum([e2[n] * grad[n] for n in range(len(qcen))], axis=0)

    if not polar_representation:        
        return tilt_comp1, tilt_comp2, e1, e2

    if polar_representation:

        # tilt angle (hope I'm not making a mistake choosing arctan2 instead of arctan)
    #     tilt_angle = np.arctan(tilt_comp2, tilt_comp1)
        tilt_angle = np.rad2deg(np.arctan2(tilt_comp2, tilt_comp1))

        # tilt magnitude
        tilt_magn = np.sqrt(tilt_comp1 **2. + tilt_comp2 **2.)

        return tilt_magn, tilt_angle, e1, e2