C-API version gvec_to_xy fails to return `nan`s in certain cases.

[donald-e-boyce]

The C-API transforms version of gvec_to_xy does not return nans under the same conditions as the numpy version.

The numpy version uses two conditions to return an nan.

1. For diffraction to occur, the g-vector needs to make an acute angle with the negative beam direction. If that condition is not met, then nans are returned.
2. If the direction of the diffracted beam points away from the detector (dot product of diffracted beam direction and detector normal is nonnegative), then again the nans are returned.

The relevant section of code is:

hexrd/hexrd/transforms/xf.py

Lines 180 to 207 in 32009f3

	# dot with beam vector (upstream direction)
	bDot = np.dot(-bHat_l.T, gVec_l).squeeze()
	# see who can diffract; initialize output array with NaNs
	canDiffract = np.atleast_1d(
	np.logical_and(bDot >= ztol, bDot <= 1.0 - ztol)
	)
	npts = sum(canDiffract)
	retval = np.nan * np.ones_like(gVec_l)
	if np.any(canDiffract):
	# subset of admissable reciprocal lattice vectors
	adm_gVec_l = gVec_l[:, canDiffract].reshape(3, npts)
	# initialize diffracted beam vector array
	dVec_l = np.empty((3, npts))
	for ipt in range(npts):
	dVec_l[:, ipt] = np.dot(
	makeBinaryRotMat(adm_gVec_l[:, ipt]), -bHat_l
	).squeeze()
	# ###############################################################
	# displacement vector calculation
	# first check for non-instersections
	denom = np.dot(nVec_l.T, dVec_l).flatten()
	dzero = abs(denom) < ztol
	denom[dzero] = 1.0 # mitigate divide-by-zero
	cantIntersect = denom > 0.0 # index to dVec_l that can't hit det

The C-API version seems to check the same conditions (although in a slightly different way), but somehow they are ignored. The relevant code is:

hexrd/hexrd/transforms/transforms_CFUNC.c

Lines 461 to 493 in 32009f3

	double bDot = -v3_v3s_dot(bHat_l, gVec_l, 1);
	double ztol = epsf;
	if ( bDot >= ztol && bDot <= 1.0-ztol ) {
	/*
	* If we are here diffraction is possible so increment the number of
	* admissable vectors
	*/
	double brMat[9];
	makeBinaryRotMat_cfunc(gVec_l, brMat);
	double dVec_l[3];
	m33_v3s_multiply(brMat, bHat_l, 1, dVec_l);
	double denom = v3_v3s_dot(nVec_l, dVec_l, 1);
	if (denom > ztol) {
	double u = num/denom;
	double v3_tmp[3];
	/* v3_tmp = P0_l + u*dVec_l - tVec_d */
	for (int j=0; j<3; j++)
	v3_tmp[j] = P0_l[j] + u*dVec_l[j] - tVec_d[j];
	result[0] = v3_v3s_dot(v3_tmp, rMat_d + 0, 3);
	result[1] = v3_v3s_dot(v3_tmp, rMat_d + 1, 3);
	/* result when computation can be finished */
	return;
	}
	}
	/* default result when computation can't be finished */
	result[0] = NAN;
	result[1] = NAN;

Note that the diffracted beam direction in the C version (dvec_l) is the negative of the diffracted beam direction in numpy (also dvec_l), but the C version checks for denom > 0 while the numpy version checks for denom < -ztol, so they are the same except for the case that the diffracted beam is (very near) parallel to the detector plane.

So I don't understand how the C-version fails to return nans in so many cases.