% This handy little code can be run on the Matlab command-line with the command: 
% showLaue(NORMALPLANE,MAXINDEX,DIST,LP,LABELALL).
% It opens a new window, which projects a picture of Laue spots as if the user
%   is the crystal looking at its reflection in the polaroid film with the incident
%   radiation wavevector opposite to the specified plane normal NORMALPLANE = [h,k,l].
% Spots are shown in all visible permutations of  [-MAXINDEX:MAXINDEX].  Projection is
%   made DIST units from crystal, DIST > 0.
% If LABELALL == 1, then all visible points are labelled with the corresponding plane.
% If LABELALL == 0, then only those visible points which belong to the form specified
%    by LP = [h,k,l] are labelled.
% Disable labelling by setting LABELALL = 0 and LP = [0,0,0].
% Axes are indicated on the image with coloured arrows: (x,y,z) = (red,green,blue).
%    So, for example, to plot spots opposite to the plane normal [1,0,1] plus or
%    minus 3, with projection distance of 12 cm, labeling only the [2,2,2] points,
%    you would type: showLaue([1,0,1],3,12,[2,2,2],0)

function showLaue(normalPlane,maxIndex,dist,LP,labelall)

%SHOWLAUE Displays a figure of Laue spots in 3-D projection

%

%  SHOWLAUE(NORMALPLANE,MAXINDEX,DIST,LP,LABELALL)

%  Projects a picture of Laue spots as if the user is the crystal

%  looking at its reflection in the polaroid film with the incident

%  radiation wavevector opposite to the specified plane normal

%  NORMALPLANE = [h,k,l].  Spots are shown in all visible permutations of

%  [-MAXINDEX:MAXINDEX].  Projection is made DIST units from crystal, DIST

%  > 0.  If LABELALL == 1, then all visible points are labelled with the

%  corresponding plane.  If LABELALL == 0, then only those visible points

%  which belong to the form specified by LP = [h,k,l] are labelled.

%  Disable labelling by setting LABELALL = 0 and LP = [0,0,0].  Axes are

%  indicated on image with coloured arrows.  (x,y,z) = (red,green,blue)

%

%  David Wong, 2004.



%  Change Log:

%  09/28/2004 - Created for utility

%  09/29/2004 - Added labelling feature

%  10/11/2004 - Added help documentation

%  



% Generate all plane G vector normals.

% Adding additional highlight flag to each vector.

temp = [0:(maxIndex*2+1)^3-1]';

Gv(:,1) = mod(floor(temp/(maxIndex*2+1)^2),maxIndex*2+1)-maxIndex;

Gv(:,2) = mod(floor(temp/(maxIndex*2+1)),maxIndex*2+1)-maxIndex;

Gv(:,3) = mod(temp,maxIndex*2+1)-maxIndex;

Gv(:,4) = zeros((maxIndex*2+1)^3,1);



% Mark the form corresponding to the label plane, LP

Gv(find((abs(Gv(:,1)) == LP(1))&(abs(Gv(:,2)) == LP(2))&(abs(Gv(:,3)) == LP(3))),4) = 1;

Gv(find((abs(Gv(:,1)) == LP(1))&(abs(Gv(:,2)) == LP(3))&(abs(Gv(:,3)) == LP(2))),4) = 1;

Gv(find((abs(Gv(:,1)) == LP(2))&(abs(Gv(:,2)) == LP(1))&(abs(Gv(:,3)) == LP(3))),4) = 1;

Gv(find((abs(Gv(:,1)) == LP(2))&(abs(Gv(:,2)) == LP(3))&(abs(Gv(:,3)) == LP(1))),4) = 1;

Gv(find((abs(Gv(:,1)) == LP(3))&(abs(Gv(:,2)) == LP(1))&(abs(Gv(:,3)) == LP(2))),4) = 1;

Gv(find((abs(Gv(:,1)) == LP(3))&(abs(Gv(:,2)) == LP(2))&(abs(Gv(:,3)) == LP(1))),4) = 1;



% Mark all of the points if labelall flag is enabled

if (labelall == 1)

    Gv(:,4) = 1;

end



% Eliminate zero vector (null!)

Gv(((maxIndex*2+1)^3+1)/2,:) = [];



% Optional: Apply structure factor condition

% BCC (sum = even)

% Fluorite (all odd, or all even)

modsum = mod(Gv(:,1)+Gv(:,2)+Gv(:,3),2);

summod = mod(Gv(:,1),2)+mod(Gv(:,2),2)+mod(Gv(:,3),2);

%missRefl = find((summod == 1)|(summod == 2));

missRefl = find(modsum == 1);

Gv(missRefl,:) = [];



ki = -normalPlane/norm(normalPlane);



% Produce set of kr vectors that represent the reflected wavevector

% directions



Gv_2 = Gv(:,1).^2 + Gv(:,2).^2 + Gv(:,3).^2;

Gv_ki = Gv(:,1:3).*ki(ones(size(Gv,1),1),:);

Gv_ki = Gv_ki(:,1) + Gv_ki(:,2) + Gv_ki(:,3);

kr(:,1) = -2*Gv(:,1) ./ Gv_2 .* Gv_ki + ki(1);

kr(:,2) = -2*Gv(:,2) ./ Gv_2 .* Gv_ki + ki(2);

kr(:,3) = -2*Gv(:,3) ./ Gv_2 .* Gv_ki + ki(3);



% Solve the locations of the wavevector lines projected onto plane governed

% by normal ki, distance dist

% Using normalized plane normal means ki dot dist = -dist.



% Calculate dot product of ki with corresponding kr

ki_kr = kr.*ki(ones(size(kr,1),1),:);

ki_kr = ki_kr(:,1) + ki_kr(:,2) + ki_kr(:,3);

badkikr = find(ki_kr == 0);

goodkikr = ki_kr;

goodkikr(badkikr) = [];

goodkr = kr;

goodkr(badkikr,:) = [];

Gv(badkikr,:) = [];

% Calculate plane unit vectors and plane magnitudes

P_mag = -dist ./ goodkikr;

P_vec = P_mag(:,ones(1,3)).*goodkr;



% Calculate G vector dot kr

Gv_kr = Gv(:,1:3).*goodkr;

Gv_kr = Gv_kr(:,1) + Gv_kr(:,2) + Gv_kr(:,3);



% Remove projection vectors for which the reflected radiation points in

% opposite direction of the plane normal vector

temp = find((Gv_kr <= 0));

badRefl = temp;

% Remove projection vectors for which the projection vector points in

% opposite direction of reflected radiation (nonsense)

temp = find(P_mag <= 0);

badRefl = union(badRefl,temp);

% Remove projection vectors lying outside of a reasonable neighbourhood of

% the origin

temp = find(P_mag > 10*dist);

badRefl = union(badRefl,temp);

P_vec(badRefl,:) = [];

Gv(badRefl,:) = [];



% Prepare labels for LP points

unlabel = find(Gv(:,4) == 0);

unlabelPoint = P_vec(unlabel,:);

label = find(Gv(:,4) == 1);

labelPoint = P_vec(label,:);

Gv = Gv(label,:);



% Find maximum of all values in set of projection vectors to set equal axis

% limits in x, y, and z.

Pmax = max(max(abs(P_vec)));



% Display 3-D image of results viewed with normal facing away user. (As if

% the user is the crystal, looking at its image projected onto the film)

[theta,phi,r] = cart2sph(-normalPlane(1),-normalPlane(2),-normalPlane(3));

theta = theta*180/pi+90;

phi = phi*180/pi;

figure;



axis equal;

grid on;

hold on;

% Show incident radiation vector

quiver3(0,0,0,ki(1),ki(2),ki(3),10,'c');

% Show x,y,z basis in r,g,b

quiver3(0,0,0,1,0,0,20,'r');

quiver3(0,0,0,0,1,0,20,'g');

quiver3(0,0,0,0,0,1,20,'b');

% Show centre of radiation (crystal)

plot3(0,0,0,'blacko');

% Plot projection points

plot3(unlabelPoint(:,1),unlabelPoint(:,2),unlabelPoint(:,3),'black.');

plot3(labelPoint(:,1),labelPoint(:,2),labelPoint(:,3),'red.');

% Label projection points

for x = 1:size(Gv,1)

    str = sprintf('%d%d%d',Gv(x,1),Gv(x,2),Gv(x,3));

    text(labelPoint(x,1),labelPoint(x,2),labelPoint(x,3),str,'FontSize',7,'HorizontalAlignment','center');

end

view(theta,phi);

axis([-Pmax,Pmax,-Pmax,Pmax,-Pmax,Pmax]);