pkposcalc_strain negative values

[MrCrab-2019]

The current function outputs negative values for lambda and ttheta for some hkl values, which are later 'filtered' out by the algorithm around line 1532. If the wavelength is less than zero, the algorithm discards it here, and builds a new array called hkl1, of only physical values.

As we are working with Q in the lab frame, real space, by the time we have applied the UB and rotation transformations, then I think it is more correct to count only those hkls with negative Qz components. Those with a positive Qz can only be satisfied where |ki| != |kf|.

In the grand scheme of things, this doesn't really matter or diamond, as the symmetry is such that there is a lot of redundancy, and it works in the end. However it is worth us checking at this point that we are discarding the correct hkls, particularly if in the future we want to correct for lower symmetry (sapphire) anvils. Malcolm, can you check if this makes sense?

I have rewritten your code a little to make it easier for me to follow (see below). If you remove the minus sign in the alp = np.degrees(np.arccos(-Q[:,2]/magQ)) line, then it reproduces your version of the code.

For convenience, the inputs are generated here:

UB = np.array([[-0.28,0.0002,0.008],[-0.007,-0.1338,-0.246],[0.0004,-0.24,-0.13]])
setang= [0,0,0]
epsilon= np.array([[1,0,0],[0,1,0],[0,0,1]])
hkl = allowedDiamRefs(-2,2,-2,2,-2,2)

Actual function:

def pkposcalc_strain(hkl, UB, setang, epsilon):
    '''
    Need to check something here regarding how negative wavelengths are dealt with later on
    Is Qz<0 then lambda is negative, and point is later disposed of (see line 1530 in old code)
    Is this correct?    

    calculates some useful numbers from (ISAW calculated) UB
    hkl is a 2D array containing all hkl's
    n.b. Adapted from the original pkposcalc to include infinitesimal strain
    present in the diamond

    The 3x3 matrix epsilon is the infinitesimal strain tensor.
    Thus, the diagonal elements are normal strains and off diagonal elements
    are shear strains.

    Here, the strain tensor operates on the reciprocal lattice vector. Accordingly
    any strain component that is > 1 implies lengthening of the corresponding
    Q component and, therefore, a real-space *compression*

    As epsilon is symmetric it includes only three independent terms
    c1,s1 - cos/sin omega
    c2,s2 - " phi
    c3,s3 - " chi
    '''

    c1, s1 = np.cos(np.radians(setang[0])), np.sin(np.radians(setang[0]))
    c2, s2 = np.cos(np.radians(setang[1])), np.sin(np.radians(setang[1]))
    c3, s3 = np.cos(np.radians(setang[2])), np.sin(np.radians(setang[2]))
    thkl = hkl.transpose()

    Rx = np.array([[1,0,0],[0,c1,-s1],[0,s1,c1]])
    Ry = np.array([[c2,0,s2],[0,1,0],[-s2,0,c2]])
    Rz = np.array([[c3,-s3,0],[s3,c3,0],[0,0,1]])
    Rall = Rz.dot(Ry).dot(Rx)

    Q = (2*np.pi*Rall.dot(epsilon.dot(UB.dot(thkl)))).transpose()
    magQ=np.empty([len(Q),])
    for i in range (0,len(Q)):
        magQ[i] = np.sqrt(Q[i].dot(Q[i]))

    d = (2*np.pi / magQ).transpose()  # by definition (note ISAW doesn't use 2pi factor)
    alp = np.degrees(np.arccos(-Q[:,2]/magQ)) #dot product of unit vector along beam and Q related to 2theta
    ttheta = (180 - 2*alp).transpose()
    lambda_1 = (2 * d * np.sin(np.radians(ttheta / 2))).transpose()
    return lambda_1, d, ttheta, Q

[mguthriem]

As we are working with Q in the lab frame, real space, by the time we have applied the UB and rotation transformations, then I think it is more correct to count only those hkls with negative Qz components. Those with a positive Qz can only be satisfied where |ki| != |kf|.

I agree with this and think it's a better way to implement the Bragg condition. I spent 30 mins going over my own code and convinced myself that what I did is essentially equivalent to this anyway (it was just needlessly complicated).

The rationale is:

lam = 2*d*sin(theta)

The relevant values of theta cannot exceed 90°, and d is always positive, therefore negative lam implies negative theta. I introduced the parameter alp. This is just the other angle in the isosceles with 2theta as the vertex and, if I cut the isosceles in half:

theta = pi/2 - alp

so, negative theta can only occur if

alp > pi/2

Since, just from geometry,

cos(alp) = Qz/Q and Q is always > 0 (in my code it's calculated directly as the product of the Q-vector with itself)

This will only be the case if Qz<0

Your approach is simpler though.

[MrCrab-2019]

Nice - ok glad that you agree. I went over it a few times myself, and came to the conclusion that it must be the case. I’ve made some good progress on putting together the other bits we discussed. Will update in due course.