%%                        MTEX Workshop 2020
%
%%                         Exercises Day 1
%
% * The exercises will be discussed during the first exercise session.
% * The exercises contain hint for useful commands. Type

help commandName

% into the command window to get additional information to this command
%
% * Exercises marked by an (*) might be a little bit more complicated.
% 
% 
%% Exercise 1 - Matlab Basics
%
%% a)
% Generate the list of squares of all numbers from 1 to 20.
%
% useful commands: .*, .^

n = (1:20).^2

%% b)
% Single out all even square numbers.
%
% useful commands: iseven, 1:2:end

n(2:2:end)

%% c)
% Compute the sum of the reciprocals of all square integers. Is there a
% limit if the numbers go to infinity?
%
% useful command: sum, 1./n

sum(1./(1:10000000).^2)- pi^2/6


%% d*)
% Compute the alternating sum of the reciprocals of the even and odd
% squares: 1/1^2 - 1/2^2 + 1/3^2 - 1/4^2 + 1/5^2 .....
% What is the limit? Hint: it is a multiple of pi^2.

n = 10000;
sum((-1).^(1:n) ./(1:n).^2) + pi^2/12


%% Exercise 2
%
%% a)
% Visualize the two vectors u1 = (1,0,2) and u2 = (2,3,-1) and determine
% the angle between these two vectors in degree.

u1 = vector3d(1,0,2)
u2 = vector3d(2,3,-1)

plot(u1)
hold on
plot(u2)
hold off

angle(u1,u2) ./ degree


%% b)
% Find a third vector u3 that is orthogonal to u1 and u2.
%
% useful command: cross

u3 = cross(u1,u2)



%% c)
% Normalize u1, u2, u3 such that they form an orthonormal basis. How can we
% check that u1, u2, u3 form a right handed basis?
%
% useful command: det

u = normalize([u1,u2,u3])

det(u1,u2,u3)


%% d)
% Find the active rotation |rot| that maps u1 to v1 and u2 to v2. Verify
% your result by computing the angle between |rot*u1| and |v1| and |rot*u2|
% and |v2|
%
% useful command: rotation.map

v1 = vector3d(1,1,1)
v2 = vector3d(1,1,-2)

rot = rotation.map(u1,v1,u2,v2)

%% e)
% Compute its rotational axis and its rotational angle of the rotation
% defined in d).


%% f)
% Could we have replaced v2 by (1,1,-1)? Why not?

v1 = vector3d(1,1,1)
v2 = vector3d(1,1,-1)

rot = rotation.map(u1,v1,u2,v2)

%% e)
% Find an approximate rotation that maps u1 onto v1 and u2 onto v2 =
% (1,1,-1). Visualize the result and determine the misfit.
%
% useful command: rotation.fit, scatter3d

rot = rotation.fit([u1,u2],[v1,v2])

angle(rot*u1,v1)./degree
angle(rot*u2,v2)./degree

scatter3d([v1,v2])
hold on
scatter3d(rot * [u1,u2])
hold off


%% Exercise 3
%
%% a)
% Define the rotation with the Euler angles (10, 20, 30) without using the
% command rotation.byEuler and check you result.
%
% useful comand: rotation.byAxisAngle


rot = rotation.byAxisAngle(zvector,10*degree) * ...
  rotation.byAxisAngle(xvector,20*degree) * ...
  rotation.byAxisAngle(zvector,30*degree)


%% b)
% Visualize this rotation in 3d Euler angle space.

plot(rot)

%% c)
% Visualize the rotation in Euler angle sections.
%
% useful commands: plotSection

plotSection(rot,'sections',12)

%% Exercise 4
%
%% a)
% Generate two random rotations rot1 and rot2 and visualize them in a three
% dimensional axis angle plot and in the three dimensional Euler angle
% plot.
% Make a nice plot with legend and filled plot markers.

rot1 = rotation.rand;
rot2 = rotation.rand;

plot(rot1,'displayName','rot1','axisAngle','filled')
hold on
plot(rot2,'displayName','rot2','axisAngle','filled')
hold off

legend show

%% b)
% Determine the misrotation angle delta between rot1 and rot2.

angle(rot1,rot2) ./ degree

%% c)
% Determine a rotation rotM such that misrotation angles between rotM and
% rot1 and rotM and rot2 are both delta/2. Highlight rotM in the two plots
% of a) as a black square.
%
% useful commands: be creative, there are many solutions :)

rotM = mean([rot1,rot2])

angle(rotM,[rot1,rot2])./degree

hold on
plot(rotM,'Marker','s','filled','MarkerFaceColor','black')
hold off

%% d*)
% Find a list rotL of 101 rotations that gradually change from rot1 to
% rot2, i.e., rotL(1) = rot1, rot(101) = rot2; rotL(51) = rotM. Visualize
% this line.
%
% useful commands: linspace

omega = linspace(0,angle(rot1,rot2),101);
misAxis = axis(rot2 * inv(rot1));
dRot = rotation.byAxisAngle(misAxis, omega);

rotL = dRot .* rot1;

hold on
line(rotL,'linewidth',2,'lineColor','green')
hold off


%% e*)
% Determine as many rotations rotP as posible such that the distances
% between (rotP,rot1) and (rotP,rot2) are equal.
%
% useful commands: vector3d/perp

pAxis = perp(misAxis);

rotP = rotation.byAxisAngle(pAxis,rand(10)*20*degree) * rotM;

hold on
plot(rotP,'markerFaceColor','orange','MarkerEdgeColor','none')
hold off

%% f*)
% Recreate the above plot in Euler angle space and Rodrigues Frank space.

%% Exercise 5
%
%% a)
% Generate random pairs of rotations with misrotation angle delta = 20
% degree.

N = 20;
rot1 = rotation.rand(N);
rot2 = rot1 .* rotation.byAxisAngle(vector3d.rand(N),20*degree);

angle(rot1,rot2) ./ degree



%% b) 
% Plot these pairs of rotations as lines in a three dimensional Euler angle
% plot-
%
% useful commands: line

for n = 1:length(rot1)
  
  line([rot1(n),rot2(n)],'line','linewidth',2)
  hold on
end
hold off

%% c) 
% Plot these pairs of rotations as lines in a three dimensional axis angle
% plot.

%% Exercise 6
%
% Consider the following Quartz crystal

cs = crystalSymmetry.load('quartz.cif')

cS = crystalShape.quartz

plot(cS,'colored')

%% a)
% Visualize the plane normals of the positive rhomboedron (1,0,-1,1), the
% negative rhomboedron (0,1,-1,1) and left tridiagonal bipyramid (1,1,-2,1)
% in a equal angle projection. Colorize them as in the crystal shape plot
% and label them nicely.
%
% useful commands: Miller
% hint: use the option upper to plot only the upper hemisphere

pr = Miller({1,0,-1,1},cs);
nr = Miller({0,1,-1,1},cs);
tb = Miller({1,1,-2,1},cs);

plot(pr,'MarkerColor',ind2color(2),'eangle','upper')
hold on
plot(nr,'MarkerColor',ind2color(3),'eangle')
plot(tb,'MarkerColor',ind2color(4),'eangle')
hold off

%% b)
% Visualize the traces of the coresponding planes.
%
% useful command: circle 

hold on
circle(pr,'lineColor',ind2color(2),'linewidth',2)
circle(nr,'lineColor',ind2color(3),'linewidth',2)
circle(tb,'lineColor',ind2color(4),'linewidth',2)
hold off


%% c)
% Determine all angles between the above planes.

angle(pr,nr) ./ degree

%% d)
% Determine the directions of the intersections between the above planes
% and visualize them in a spherical projection.

hold on
plot(cross(pr,nr),'antipodal','MarkerColor','black')
plot(cross(nr,tb),'antipodal','MarkerColor','black')
plot(cross(tb,pr),'antipodal','MarkerColor','black')
hold off

%% Exercise 7
% 
%% a) 
% Generate a symmetry group that consist of three fold axis in direction
% (111) and a two fold axis perpendicular to it. Visualize the symmetry
% group using the command plot.
%
% useful commands: rotation.byAxisAngle, crystalSymmetry.byElements

rot3 = rotation.byAxisAngle(vector3d(1,1,1),120*degree);
rot2 = rotation.byAxisAngle(vector3d(1,-1,0),180*degree);

cs = crystalSymmetry.byElements([rot3,rot2])

plot(cs)

%% b)
% Turn the above symmetry into a Laue group by adding the inversion and
% plot it again.
%
% useful commands: crystalSymmetry/add, rotation.inversion

cs.add(rotation.inversion)

plot(cs.Laue,'lower')

%% c)
% Add to the above symmetry a four fold axis in direction (001). Which
% symmetry group do we obtain?

cs = cs.add(rotation.byAxisAngle(vector3d.Z,90*degree))

plot(cs)

%% d)
% Verify, that a mirroring at a plane with normal (111) is nothing else
% then a rotation around the axis (111) about 180 degree followed by an
% inversion.
%
% useful commands: reflection

reflection(vector3d(1,1,1))

- rotation.byAxisAngle(vector3d(1,1,1),180*degree)

%% Exercise 8
%
%% a) 
% Define trigonal crystal symmetry 321 and a random crystal direction.
% Compute all symmetrically equivalent crystal directions without using the
% command symmetrise and check that your result is correct.
%
% useful commands: Miller.rand


cs = crystalSymmetry('321')

m = Miller.rand(cs)

mSym = m.symmetrise
cs * m


%% b)
% Compute the angles between all symmetrically equivalent directions
%
% useful commands: angle, angle_outer, symmetrise
% useful options: noSymmetry

delta = angle_outer(mSym,mSym,'noSymmetry')./degree

%% c) 
% Find a crystal direction which has exactly two symmetrically equivalent
% directions with respect to the symmetry 321
%
% useful commands: unique
% useful options: noSymmetry

m = Miller(0,0,1,cs)
unique(m.symmetrise,'noSymmetry')


%% d) 
% Find a crystal direction which has exactly three symmetrically equivalent
% directions with respect to the symmetry 321

m = Miller(2,-1,-1,0,cs,'uvw')
unique(m.symmetrise,'noSymmetry')


%% e)
% What are the common symmetries of the point groups 321 and m-3?
% Do the common symmetries again form a point group?
%
% useful commands: symmetry/disjoint

cs2 = crystalSymmetry('m-3')

plot(disjoint(cs2,cs))

%% Exercises 9
%
% Define a crystal symmetry such that (011) is not parallel to {011}.
% Verify this using the command angle. What is the highest symmetry you
% could have used?

cs = crystalSymmetry('622',[1 1 3])
angle(Miller(0,1,1,cs,'hkl'),Miller(0,1,1,cs,'uvw'))

