Laue diffraction matlab code

Back to scattering techniques animations

Download matlab code

Laue2D.m

Download matlab code

Laue3D.m

% Program to generate 2D schematic movie of Laue method

clear; close all;

vid = VideoWriter('Laue2D.mp4','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);

kV1 = 1/1; % Radius of Ewald sphere, here 1 AA^-1

kV2 = 1/4; % Radius of Ewald sphere, here 1 AA^-1

% Create a set of diffraction maxima with separations of 1/(nBP) AA^-1 out to

% 1 AA (1 AA^-1)

bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))

rBP = 0.0053; % Radius of Bragg peak as plotted on ES

numBP = 14; % Number of Bragg peaks to edge of diffraction pattern at 0.5 AA^-1

% Vertices of stretched octahedral crystal

rhSize = 0.025;

v1 = [rhSize,0,0];

v2 = [0,rhSize,0];

v3 = [0,0,rhSize*1.25];

v4 = [-rhSize,0,0];

v5 = [0,-rhSize,0];

v6 = [0,0,-rhSize*1.25];

[x,y,z] = sphere; % Create 3D array to plot spheres

[a,b,c] = sphere; % Create 3D array to plot spheres

for phi = 0:0.5:179.875 % Rotate Bragg peaks around y-axis in steps of 0.5 degrees

hold off

% Plot semitransparent Ewald sphere at (-kV1, 0, 0)

[x,y,z] = sphere(70); surf(kV1*x-kV1,kV1*y,kV1*z,'FaceAlpha',0.07,'FaceColor',...

bpCol,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

hold on

% Plot semitransparent Ewald sphere at (-kV2, 0, 0)

[x,y,z] = sphere(70); surf(kV2*x-kV2,kV2*y,kV2*z,'FaceAlpha',0.08,'FaceColor',...

'r','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

% Draw crystal and rotate it

f1 = fill3(v1,v2,v3,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

f2 = fill3(v4,v2,v3,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

f3 = fill3(v1,v2,v6,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

f4 = fill3(v4,v2,v6,'green','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

rotate(f1,[0 1 0],phi,[0,0,0]);

rotate(f2,[0 1 0],phi,[0,0,0]);

rotate(f3,[0 1 0],phi,[0,0,0]);

rotate(f4,[0 1 0],phi,[0,0,0]);

% Loop through hkl space

for h = -1:1/numBP:1

for l = -1:1/numBP:1

hl = (h^2 + l^2)^0.5;

if (hl <= 1) % Only plot out Bragg peaks in reciprocal lattice

% within radius of 1 AA^-1.

hold on

LL = 1.01; % Defines limits of plotted figure

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-LL LL]);

% Rotate reciprocal lattice to angle phi

% and calculate new absolute (h k l) after rotation

hnew = h*cosd(phi) + l*sind(phi);

lnew = l*cosd(phi) - h*sind(phi);

innerIneq = (kV2 + hnew)^2 + lnew^2;

outerIneq = (kV1 + hnew)^2 + lnew^2;

if (innerIneq < kV2^2) ||... % Inside small Ewald sphere

(outerIneq > kV1^2) ||... % Outside large Ewald sphere

(hnew > 0) % All detectable BPs have h <= 0

% Undetected Bragg peaks outside Ewald volume

H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(1.0 - hl)/1.0,...

'FaceColor',bpCol,'LineStyle','none');

rotate(H,[0 1 0],phi,[0,0,0]);

elseif (hl == 0) % (000) direct beam

H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(1.0 - hl)/1.0,...

'FaceColor',bpCol,'LineStyle','none');

rotate(H,[0 1 0],phi,[0,0,0]);

else % Detected Bragg peaks lying in Ewald volume

H = surf(2*rBP*x+h,2*rBP*y,2*rBP*z+l,'FaceAlpha',...

((1.0 - (hl)^2)/1.0),...

'FaceColor','r','LineStyle','none','FaceLighting','Flat');

rotate(H,[0 1 0],phi,[0,0,0]);

end % End of finding out if BP on ES

end % End of selecting only hl values within radius of 0.5 r.l.u.

end % End of l-loop

end % End of h-loop

% End of drawing diffraction maxima

% Cross at ES origin and dotted axis

line1 = plot3([-1 -1],[-0.05 0.05],[0 0],'color','red', 'LineWidth', 1.4);

line2 = plot3([-1 -1],[0 0],[-0.05 0.05],'color','red', 'LineWidth', 1.4);

line3 = plot3([-1.05 1.05],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');

set(gca, 'Projection','perspective');

set(gca,'View',[0,0]);

lp = [-0.7 -0.5 0.5];

light('Position',lp,'Style','infinite');

axis equal

axis off

LL = 1.01;

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-LL LL]);

% Store the frame

frame = getframe(gcf);

writeVideo(vid,frame);

end % End of phi-rotation loop

% Output the movie as an mpg file

close(vid);

% Program to generate 3D schematic movie of Laue diffraction while rotating

% the crystal. The "diffracting volume" or "Ewald volume" is defined by two

% limits to the k-vector. This volume is hence bounded by a large Ewald

% sphere of radius kV1 and within this a small Ewald sphere with radius

% kV2; they touch one another at (0,0,0). As default, these are set to 1

% and 0.25 reciprocal Angstroms (12.4 keV and 3.1 keV, respectively).

% Parameter list

% ______________

% kV1 = k-vector of outer surface of Ewald volume in units of 2pi reciprocal Angstroms

% kV2 = k-vector of inner surface of Ewald volume in units of 2pi reciprocal Angstroms

% rBP = radius of spheres representing Bragg peaks in units of 2pi reciprocal Angstroms

% numBP = +/- range of BPs (max # peaks in a row = 2*numBP+1)

% bpExtentRad = extent (radius) of sphere to which Bragg peaks are plotted. Also in units of 2pi reciprocal Angstroms

% phi = rotation angle around y-axis in degrees

% h, k, l = orthogonal reciprocal-space coordinates

% hnew, knew, lnew = angularly transformed positions of (h, k, l)

% thnow = absolute Bragg angle associated with current (hnew, knew, lnew) Bragg peak

% kVnow = k-vector associated with current (hnew, knew, lnew) Bragg peak in units of 2pi reciprocal Angstroms

clear; close all;

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

vid = VideoWriter('Laue.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 30;

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

% Create semitransparent Ewald sphere centered around (-kV,0,0) with 1

% reciprocal Angstrom radius

kV1 = 1/1; % Radius of outer surface of Ewald volume, here 1 AA^-1

kV2 = 1/2.5; % Radius of inner surface of Ewald volume, here 0.4 AA^-1

% Create a set of diffraction maxima out to bpExtentRad reciprocal Angstroms

bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))

lyel = [1 1 0.64]; % Color of Bragg peaks on detector

rBP = 0.005; % Radius of Bragg peak

numBP = 10; % From -numBP to + numBP within range of +/- bpExtentRad

bpExtentRad = 0.4; % Radial extent to which Bragg peaks are shown

stepSize = bpExtentRad/numBP; % Separation of adjacent Bragg peaks in reciprocal space

% Vertices of crystal drawn around (0,0,0) point

rhSize = 0.05;

X1 = [0 0 rhSize]; % x-coordinates of vertices of first face of octahedron

Y1 = [0,-rhSize,0]; % y-coordinates of vertices of first face of octahedron

Z1 = [rhSize,0 0]; % z-coordinates of vertices of first face of octahedron

X2 = [0 0 rhSize]; % etc, etc...

Y2 = [0,rhSize,0];

Z2 = [rhSize,0 0];

X3 = [0 0 -rhSize];

Y3 = [0,-rhSize,0];

Z3 = [rhSize,0 0];

X4 = [0 0 -rhSize];

Y4 = [0,rhSize,0];

Z4 = [rhSize,0 0];

X5 = [0 0 rhSize];

Y5 = [0,-rhSize,0];

Z5 = [-rhSize,0 0];

X6 = [0 0 rhSize];

Y6 = [0,rhSize,0];

Z6 = [-rhSize,0 0];

X7 = [0 0 -rhSize];

Y7 = [0,-rhSize,0];

Z7 = [-rhSize,0 0];

X8 = [0 0 -rhSize];

Y8 = [0,rhSize,0];

Z8 = [-rhSize,0 0];

detPos = 0.8; % Position of detector face in reciprocal space

[x,y,z] = sphere; % Used for Ewald sphere

[a,b,c] = sphere; % Used for BPs on detector

[Z,Y,X] = cylinder(1,180); % Side walls of detector

detector = zeros(3,101);

th = 0:pi/50:2*pi;

detector(1,:) = zeros(1,101) + detPos; % x-coordinates of detector body

detector(2,:) = 1*cos(th); % y-coordinates

detector(3,:) = 1*sin(th); % z-coordinates

phiStep = 0.5;

for phi = 0:phiStep:89.999 % Rotate Bragg peaks around y-axis in phiStep-deg steps

phi % Monitor progress of program

laueInt = NaN((2*numBP+1)^3,5); % Array to plot as scatter3 the BPs within the Ewald volume

undetInt = NaN((2*numBP+1)^3,5); % Array to plot as scatter3 the BPs outside the Ewald volume

hold off

% Draw outer spherical surface of Ewald volume centered at x = -kV1 and

% radius kV1

[x,y,z] = sphere(70); surf(kV1*x-kV1,kV1*y,kV1*z,'FaceAlpha',0.125,'FaceColor',...

'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

hold on

% Draw inner spherical surface of Ewald volume centered at x = -kV2 and

% radius kV2

[x,y,z] = sphere(70); surf(kV2*x-kV2,kV2*y,kV2*z,'FaceAlpha',0.125,'FaceColor',...

'r','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

% Draw octahedral crystal

f1 = fill3(X1,Y1,Z1,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f2 = fill3(X2,Y2,Z2,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f3 = fill3(X3,Y3,Z3,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f4 = fill3(X4,Y4,Z4,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f5 = fill3(X5,Y5,Z5,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f6 = fill3(X6,Y6,Z6,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f7 = fill3(X7,Y7,Z7,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

f8 = fill3(X8,Y8,Z8,'green','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

rotate(f1,[0 1 0],phi,[0,0,0]);

rotate(f2,[0 1 0],phi,[0,0,0]);

rotate(f3,[0 1 0],phi,[0,0,0]);

rotate(f4,[0 1 0],phi,[0,0,0]);

rotate(f5,[0 1 0],phi,[0,0,0]);

rotate(f6,[0 1 0],phi,[0,0,0]);

rotate(f7,[0 1 0],phi,[0,0,0]);

rotate(f8,[0 1 0],phi,[0,0,0]);

hklstep = (-bpExtentRad:stepSize:bpExtentRad);

ii = 0;

for h = hklstep

for k = hklstep

for l = hklstep

ii = ii+1;

hkl = (h^2 + k^2 + l^2)^0.5;

if (hkl <= bpExtentRad) % Plot out Bragg peaks within radius of bpExtentRad r.l.u.

hold on

LL = 1.05;

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-LL LL]);

% Calculate new absolute h k l after rotation

hnew = h*cosd(phi) + l*sind(phi);

knew = k;

lnew = l*cosd(phi) - h*sind(phi);

% innerIneq and outerIneq define limits of BPs that are

% detected

innerIneq = (kV2 + hnew)^2 + knew^2 + lnew^2;

outerIneq = (kV1 + hnew)^2 + knew^2 + lnew^2;

if (innerIneq < kV2^2) ||... % Inside small Ewald sphere

(outerIneq > kV1^2) ||... % Outside large Ewald sphere

(hnew > 0) % All detectable BPs have h <= 0

% Undetected Bragg peaks outside Ewald volume

undetInt(ii,1) = hnew;

undetInt(ii,2) = knew;

undetInt(ii,3) = lnew;

else

% Detected Bragg peaks inside Ewald volume

laueInt(ii,1) = hnew;

laueInt(ii,2) = knew;

laueInt(ii,3) = lnew;

% Determine size of present k-vector

thnow = atan2(abs(hnew),(knew^2 + lnew^2)^0.5); % theta for this BP

kVnow = abs(hkl/(2*sin(thnow))); % Present k-vector for this BP

% Plot Bragg peaks detected on detector

apos = detPos - 0.01;

bpos = knew * (kVnow+detPos)/(kVnow-abs(hnew));

cpos = lnew * (kVnow+detPos)/(kVnow-abs(hnew));

BPpos = (bpos^2 + cpos^2)^0.5; % Radius of position of BP on detector

if (BPpos < 1.0) % Plot BPs on detector face

detBPsize = 2*((kV1 - kVnow)/kV1)*rBP;

BP = surf(0.02*rBP*a+apos,detBPsize*b+bpos,detBPsize*c+cpos,...

'FaceAlpha',1, ...

'FaceColor',lyel,'LineStyle','none');

end % End of plotting BP signal on detector

end % End of finding out if BP in Ewald Volume

end % End of selecting only hkl values within radius of bpExtentRad

end % End of l-loop

end % End of k-loop

end % End of h-loop

% End of drawing diffraction maxima

scatter3(laueInt(:,1),laueInt(:,2),laueInt(:,3),3.0,'r',...

'MarkerFaceColor','Flat','MarkerEdgeColor','r');

scatter3(undetInt(:,1),undetInt(:,2),undetInt(:,3),3.0,'b',...

'MarkerFaceColor','Flat','MarkerEdgeColor','b');

% Side wall of detector

surf(detPos+X/5,Y*1.05,Z*1.05,'FaceAlpha',1,'FaceColor',[0.1 0.1 0.1],...

'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

% Detector face

fill3(detector(1,:),1.05*detector(2,:),1.05*detector(3,:),[0 0 0]);

fill3(detector(1,:)-0.005,detector(2,:),detector(3,:),[0.2 0.1 0.1]);

% Side wall of beamstop

surf(detPos-0.1+X/12.5,0.03*Y,0.03*Z,'FaceAlpha',1,'FaceColor',[0.2 0.2 0.2],...

'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

% Beamstop face

fill3(detector(1,:)-0.1,0.03*detector(2,:),0.03*detector(3,:),[0.4 0.4 0.4],...

'LineStyle','none');

% Cross at outer ES origin and dotted axis

line1 = plot3([-1 -1],[-0.05 0.05],[0 0],'color','red', 'LineWidth', 1.4);

line2 = plot3([-1 -1],[0 0],[-0.05 0.05],'color','red', 'LineWidth', 1.4);

line3 = plot3([-1.05 detPos-0.1],[0 0],[0 0],'color','red', 'LineWidth', 1.4, 'LineStyle','-.');

set(gca, 'Projection','perspective');

lp = [-0.7 -0.5 0.5];

light('Position',lp,'Style','infinite');

set(gca,'View',[-40,20]);

axis equal

axis off

LL = 1.05*kV1;

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-LL LL]);

% Store the frame

frame = getframe(gcf);

writeVideo(vid,frame);

end % End of phi-rotation loop

% Output the movie as an mpg file

close(vid);

% Program to generate 2D schematic movie of Laue method with detectable

% Bragg peaks color coded according to the photon energy to which they

% correspond. After the first 180-degree rotation, a high-pass energy

% filter is applied in the detector that varies its cut-off energy from that

% of the inner Ewald sphere to that of the outer sphere over the next 180

% degrees, causing the Bragg peaks within the Ewald volume to slowly

% vanish as the cutoff energy exceeds their photon energy. This is an

% important feature of modern detectors, enabling the "overlap problem"

% associated with Laue diffraction to be overcome.

%

clear; close all;

vid = VideoWriter('Laue2DcolorCoded4energy.mp4','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);

kV1 = 1/1; % Radius of Ewald sphere, here 1 AA^-1

kV2 = 1/4; % Radius of Ewald sphere, here 1 AA^-1

% Create a set of diffraction maxima with separations of 1/(nBP) AA^-1 out to

% 1 AA (1 AA^-1)

bpCol = [0.3 0.28 0.55]; % Color of Bragg peaks (if not on surface of Ewald sphere (ES))

rBP = 0.0053; % Radius of Bragg peak as plotted on ES

numBP = 14; % Number of Bragg peaks to edge of diffraction pattern at 0.5 AA^-1

rhSize = 0.04;

% red - yellow - green - blue - violet

rgb = customcolormap([0 0.25 0.5 0.75 1], {'#ff0000','#ffff00','#00ff00','#0000ff','#7f00ff'});

fliprgb = flip(rgb);

colormap(fliprgb);

[x,y,z] = sphere; % Create 3D array to plot spheres

[a,b,c] = sphere; % Create 3D array to plot spheres

for phi = 0:0.5:359.875 % Rotate Bragg peaks around y-axis in steps of 0.5 degrees

hold off

% Plot semitransparent Ewald sphere at (-kV1, 0, 0)

[x,y,z] = sphere(70); surf(kV1*x-kV1,kV1*y,kV1*z,'FaceAlpha',0.07,'FaceColor',...

[0.5 0 1],'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

hold on

% Plot semitransparent Ewald sphere at (-kV2, 0, 0)

[x,y,z] = sphere(70); surf(kV2*x-kV2,kV2*y,kV2*z,'FaceAlpha',0.08,'FaceColor',...

'r','LineStyle','none','FaceLighting','flat','DiffuseStrength',1);

% Draw crystal and rotate it

[V,F] = platonic_solid(3,rhSize); % Octahedron, size = 1

ps = patch('Faces',F,'Vertices',V,'FaceColor','g','FaceAlpha',1, ...

'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1);

direction = [0 1 0];

rotate(ps,direction,phi,[0 0 0])

% Loop through hkl space

for h = -1:1/numBP:1

for l = -1:1/numBP:1

hl = (h^2 + l^2)^0.5;

if (hl <= 1) % Only plot out Bragg peaks in reciprocal lattice

% within radius of 1 AA^-1.

hold on

LL = 1.01; % Defines limits of plotted figure

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-LL LL]);

% Rotate reciprocal lattice to angle phi

% and calculate new absolute (h k l) after rotation

hnew = h*cosd(phi) + l*sind(phi);

lnew = l*cosd(phi) - h*sind(phi);

%hl = (hnew^2 + lnew^2)^0.5;

innerIneq = (kV2 + hnew)^2 + lnew^2;

outerIneq = (kV1 + hnew)^2 + lnew^2;

if (innerIneq < kV2^2) ||... % Inside small Ewald sphere

(outerIneq > kV1^2) ||... % Outside large Ewald sphere

(hnew > 0) % All detectable BPs have h <= 0

% Undetected Bragg peaks outside Ewald volume

H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(1.0 - hl)/1.0,...

'FaceColor',bpCol,'LineStyle','none');

rotate(H,[0 1 0],phi,[0,0,0]);

elseif (hl == 0) % (000) direct beam

H = surf(rBP*x+h,rBP*y,rBP*z+l,'FaceAlpha',(1.0 - hl)/1.0,...

'FaceColor',bpCol,'LineStyle','none');

rotate(H,[0 1 0],phi,[0,0,0]);

else % Detected Bragg peaks lying in Ewald volume

% Determine size of present k-vector

thnow = atan2(abs(hnew),abs(lnew)); % theta for this BP

kVnow = abs(hl/(2*sin(thnow))); % Present k-vector for this BP

colorIndex = round(256*(kVnow-kV2)/(kV1-kV2))+1;

if (colorIndex==257)

colorIndex = 256;

end

if (phi>=180) % Apply high-pass filter that varies from

% energy corresponding to kV2 to that

% corresponding to kV1 over the next 180

% degrees

kVcutoff = (kV2 + (phi-180)*(kV1-kV2)/180);

if (kVnow<=kVcutoff)

BPalpha = 0; % BP disappears!

else

BPalpha = ((1.0 - (hl)^2)/1.0);

end

else

BPalpha = ((1.0 - (hl)^2)/1.0);

kVcutoff = kV2;

end

H = surf(2*rBP*x+h,2*rBP*y,2*rBP*z+l,'FaceAlpha',...

BPalpha,...

'FaceColor',fliprgb(colorIndex,:),'LineStyle','none','FaceLighting','Flat');

rotate(H,[0 1 0],phi,[0,0,0]);

end % End of finding out if BP on ES

end % End of selecting only hl values within radius of 0.5 r.l.u.

end % End of l-loop

end % End of h-loop

% End of drawing diffraction maxima

% Cross at ES origin and dotted axis

line1 = plot3([-kV2 -kV2],[0 0],[-0.025 0.025],'color','red', 'LineWidth', 1.76);

line2 = plot3([-kV1 -kV1],[0 0],[-0.025 0.025],'color','red', 'LineWidth', 1.76);

line3 = plot3([-kV1-0.025 kV1+0.025],[0 0],[0 0],'color','red', 'LineWidth', 1.76, 'LineStyle','-.');

% Draw dot-dashed circle corresponding to Ewald sphere at cut-off

% energy

kVcutoffPhi = 0:pi/90:2*pi;

if (phi>=180)

%kVcutoff = (kV2 + (phi-180)*(kV1-kV2)/180);

colorIndex2 = round(256*(kVcutoff-kV2)/(kV1-kV2))+1;

if (colorIndex2==257)

colorIndex2 = 256;

end

kVcutoffCoords(1,:) = kVcutoff*sin(kVcutoffPhi)-kVcutoff;

kVcutoffCoords(2,:) = 0.0*sin(kVcutoffPhi);

kVcutoffCoords(3,:) = kVcutoff*cos(kVcutoffPhi);

plot3(kVcutoffCoords(1,:),kVcutoffCoords(2,:),kVcutoffCoords(3,:), ...

'color',fliprgb(colorIndex2,:), 'LineWidth', 1.76, 'LineStyle','-.');

end

cbar = colorbar;

set(cbar,'position',[.2 .23 .01 .2]); % Size and position of colorbar

cbar.Color = [0 0 0];

caxis([kV2 kV1])

set(cbar,'Ticks',[kV2,kV1],...

'TickLabels',{'h\nu_1','h\nu_2'},'FontSize',16);

set(gca, 'Projection','perspective');

set(gca,'View',[0,0]);

lp = [-0.7 -0.5 0.5];

light('Position',lp,'Style','infinite');

axis equal

axis off

LL = 1.01;

xlim([-LL LL]);

ylim([-LL LL]);

zlim([-LL LL]);

% Store the frame

frame = getframe(gcf);

writeVideo(vid,frame);

end % End of phi-rotation loop

% Output the movie as an mpg file

close(vid);