Generation of lunes on the Ewald sphere
Matlab code
% Program to generate lunes
% Uses an adaptation of the function platonic_solid, see
% https://ch.mathworks.com/matlabcentral/fileexchange/28213-platonic_solid
% Courtesy Kevin Moerman
clear; close all;
prompt = 'Rotation of crystal around z-axis in degrees [default = 0]: ';
theta = input(prompt);
if isempty(theta)
theta = 0;
end
fileName = 'lunes.mp4'; % Name of file depends on type of crystal sample
vid = VideoWriter(fileName,'MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
set(gca, 'Projection','perspective');
lp = [-1 -2 0.7];
myBlue = [0.55 0.61 0.89];
myBlue2 = [0.31 0.37 0.64];
myGold = [0.7969 0.6406 0.2383];
LL = 1.1;
delk = 0.01;
[a,b,c] = sphere(70);
for phi = 0:0.25:89.8
phi
newplot
hold off
xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);
axis equal
axis square
axis off
for hh = -1:0.02:1
for kk = -1:0.02:1
for ll = -1:0.02:1
% Calculate new absolute h k l after rotation
hhnew = hh*cosd(phi) + ll*sind(phi);
kknew = kk;
llnew = ll*cosd(phi) - hh*sind(phi);
hhnew2 = hhnew*cosd(theta) + kknew*sind(theta);
kknew2 = kknew*cosd(theta) - hhnew*sind(theta);
llnew2 = llnew;
radNow = ((hhnew2-0.5)^2+kknew2^2+llnew2^2)^0.5;
if (radNow <= 0.5+delk/2.5) && (radNow >= 0.5-delk/2.5)
[V,F] = platonic_solid(2,delk/2); % Cube
V(:,1) = V(:,1)+hhnew2;
V(:,2) = V(:,2)+kknew2;
V(:,3) = V(:,3)+llnew2;
ps = patch('Faces',F,'Vertices',V,'FaceColor',myGold,'FaceAlpha',1, ...
'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1);
direction = [0 1 0];
rotate(ps,direction,phi,[hhnew2 kknew2 llnew2])
hold on
end
end
end
end
xlim([-0.1 LL]);
ylim([-LL LL]/1.1);
zlim([-LL LL]/1.1);
% Ewald sphere
surf(0.5*a+0.5,0.5*b,0.5*c,'FaceAlpha',0.55,'FaceColor',...
myBlue2,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.25);
% Dummy sphere to stop small changes in FOV
surf(0.52*a+0.5,0.52*b,0.52*c,'FaceAlpha',0.0001,'FaceColor',...
myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
axis equal
axis off
axis tight
set(gca, 'Projection','perspective');
set(gca,'View',[-30,32]);
light('Position',lp,'Style','infinite');
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end
% Output the movie as an mpg file
close(vid);