gtCrystHKL2CartesianMatrix: made the function more simple to read and understand

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Cryst/gtCrystHKL2CartesianMatrix.m

@@ -12,8 +12,10 @@ function [Bmat, Amat] = gtCrystHKL2CartesianMatrix(lp)
 %       - Y lies in the a1-a2 plane
 %       - Z is cross(X,Y)
 % 
-%     The definition corresponds to the one in Poulsen's 3DRXRD book, Chapter
-%     3. This function essentially computes the B matrix in the equation 
+%     The definition does not correspond to the one in Poulsen's 3DXRD
+%     book, Chapter 3, but it is an equivalent representation based on the
+%     real space lattice.
+%     This function essentially computes the B matrix in the equation 
 %     Gc = B*Ghkl  (used with column vectors!)
 %
 %     INPUT:
@@ -22,48 +24,48 @@ function [Bmat, Amat] = gtCrystHKL2CartesianMatrix(lp)
 %
 %     OUTPUT:
 %       Bmat = the transformation matrix for column vectors (!)
-%       Amat = the A matrix too (real space) for dealing with UVW 
-%              lattice directions.
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Lattice vectors in real space, in a Cartesian basis
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Calculate 'a1' and 'a2' lattice vector X,Y,Z coordinates:
-a1 = [lp(1)             0                 0];
-a2 = [lp(2)*cosd(lp(6)) lp(2)*sind(lp(6)) 0];
-
-% Calculate 'a3' components in X-Y plane:
-x3 = lp(3)*cosd(lp(5));
-y3 = lp(3)*(cosd(lp(4))-cosd(lp(5))*cosd(lp(6)))/sind(lp(6));
-
-% Calculate 'a3' Z coordinate from length
-z3 = sqrt(lp(3)^2 - x3^2 - y3^2);
+%       Amat = the A matrix too (real space) for dealing with UVW lattice
+%              directions.
+%
 
-% Assemble 'a3'
-a3 = [x3 y3 z3];
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%% Lattice vectors in real space, in a Cartesian basis
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    cos_alpha = cosd(lp(4));
+    cos_beta  = cosd(lp(5));
+    cos_gamma = cosd(lp(6));
+    sin_gamma = sind(lp(6));
 
+    % Calculate 'a', and 'b' lattice vectors in X, Y, Z coordinates:
+    a = lp(1) * [1,         0,         0];
+    b = lp(2) * [cos_gamma, sin_gamma, 0];
 
-% Volume of spanned by 'a'-s (determinant of matrix of a-s)
-crossa2a3 = [a2(2)*a3(3)-a2(3)*a3(2); a2(3)*a3(1)-a2(1)*a3(3); a2(1)*a3(2)-a2(2)*a3(1)];
-avol = a1*crossa2a3;
+    % Calculate 'c1' and 'c2' components:
+    c1 = lp(3) * cos_beta;
+    c2 = lp(3) * (cos_alpha - cos_gamma * cos_beta) / sin_gamma;
+    % Calculate 'c3' component as the result:
+    c = [c1,  c2, sqrt(([lp(3), c1, c2] .^ 2) * [1; -1; -1])];
 
-% for output
-Amat=[a1' a2' a3'];
+    % We keep it transposed for the computation of B, and then transpose it
+    % back again in the end.
+    Amat = [a; b; c];
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Reciprocal lattice vectors in the Cartesian basis
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Reciprocal lattice vectors by definition:
-%  b1 = cross(a2,a3)/avol
-%  b2 = cross(a3,a1)/avol
-%  b3 = cross(a1,a2)/avol
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%% Reciprocal lattice vectors in the Cartesian basis
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Reciprocal lattice vectors by Crystallographers' definition:
+    %  a* = cross(b, c) / cell_vol
+    %  b* = cross(c, a) / cell_vol
+    %  c* = cross(a, b) / cell_vol
+    cross_prods = gtMathsCross(Amat, Amat([2 3 1], :))';
 
-b1 = crossa2a3/avol;
-b2 = [a3(2)*a1(3)-a3(3)*a1(2); a3(3)*a1(1)-a3(1)*a1(3); a3(1)*a1(2)-a3(2)*a1(1)]/avol;
-b3 = [a1(2)*a2(3)-a1(3)*a2(2); a1(3)*a2(1)-a1(1)*a2(3); a1(1)*a2(2)-a1(2)*a2(1)]/avol;
+    % Volume of the direct lattice unit cell: V = (a dot (b cross c))
+    cell_vol = a * cross_prods(:, 2);
 
-Bmat = [b1 b2 b3];
+    % B matrix finally
+    Bmat = cross_prods(:, [2 3 1]) / cell_vol;
 
+    % A matrix (for output)
+    Amat = Amat';
 
 end % of function