Powder-diffraction experiment matlab code

Back to scattering techniques animations

Download matlab code

% Cartoon movie of a modern powder diffraction experiment employing a

% microstrip detector for parallel data acquisition. The generated powder

% pattern is of silicon at 1-Angstrom radiation

clear; close all;

% -----------------------------------------------------------------------

% Storyboard of the movie

% _______________________

% 1) Let sample holder accelerate to full speed

% 2) Fade in incident radiation and Debye-Scherrer cones/intersections

% 3) Continue for a short while

% 4) Fade out full cones, leaving only the parts of the cones incident on

% the detector, plus the bright lines on detector where the rings intersect

% 5) Fade in curved diffraction pattern in front of detector and hold for a

% short while

% 6) Straighten diffraction pattern to lie on x-axis

% 7) Show Cartesian plot for a short while

% -----------------------------------------------------------------------

% Timing of movie

rampSpeed = 0.4; % degrees increase/frame during sample angular acceleration

topSpeed = 24; % top speed of rotation of sample holder/frame

sampleAccelFrames = round(topSpeed/rampSpeed) + 1; % Frames to get capillary to full speed

accel = 0:rampSpeed:topSpeed; % Values of steps in position during sample acceleration

fadeInRadiation = 60; % Number of frames to fade in the incident radiation and DS-cones

initialSitch = 60; % Number of frames to accustomize viewer to the full scattering scenario

fadeOutFullFadeInDet = 60; % Number of frames to fade in detector signal and fade out full DS-cones

fadeInRadial = 60; % Number of frames to fade in radial diffraction pattern

keepRadial = 60; % Number of frames to hold on radial diffraction pattern

radToLin = 60; % Number of frames to straighten diffraction pattern to being linear

linDiff = 90; % Number of frames showing linear diffraction pattern

constSpeedSteps = fadeInRadiation + initialSitch + fadeOutFullFadeInDet + ...

fadeInRadial + keepRadial + radToLin + linDiff; % Number of frames at top sample speed

totSteps = constSpeedSteps + sampleAccelFrames; % Total number of frames

constSpeed = topSpeed*ones(1,constSpeedSteps); % Vector of angle changes for each frame when running at full speed

totAngleSteps = [accel constSpeed]; % Vector of angle changes for each and every frame

% Form factor parameters for tetravalent silicon

a1 = 5.66269;

b1 = 2.6652;

a2 = 3.07164;

b2 = 38.6634;

a3 = 2.62446;

b3 = 0.916946;

a4 = 1.3932;

b4 = 93.5458;

c = 1.24707;

% Set up and initialize x-y data for plot out of powder intensity graph

mythenMax = 127; % Maximum 2theta angle measured by microstrip detector

tthStep = 0.025; % Step size in 2theta in degrees

powderPattx = (mythenMax-120:tthStep:mythenMax); % Range of 2theta values

powderPatty = 0.0*powderPattx; % Initialize diffraction intensities to zero

sigma = 0.16; % Standard deviation of diffraction-peak widths

coneRad = 7; % Radius of Debye-Scherrer cone

aSi = 5.431; % u.c. size of cubic silicon in Angstroms

hv = 12.3984; % Photon energy for 1 Angstrom radiation

lambda = 1.0; % X-ray wavelength

% -------------------------------------------------------

% Calculate Bragg-peak intensities and angles for silicon

hklmax = 10; % Maximum Miller index

nBPs = 0; % Initialize the number of Bragg peaks to zero

% Now work out SFs for 1/32 of volume: only positive h, k, l; k <= h, l <= k

for h = 1:hklmax

for k = 0:h

for l = 0:k

% dhkl = 1/Y. Determine here Y

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

hkl2 = (h^2 + k^2 + l^2);

Y = (hkl)/aSi; % 1/d_hkl

twoth = 2*asind(lambda*Y/2); % 2theta in degrees

% Only allowed reflections for diamond structure

% For h, k, l all odd, or h, k, l all even and h+k+l = 4m

% See plot of allowed reflections v (h^2 + k^2 + l^2) in

% Figure 6.13 of the second edition of "Introduction to

% Synchrotron Radiation - Techniques and Applications" by

% Philip Willmott, Wiley (2019)

if (mod(hkl2,8) == 0) || (mod(hkl2,8) == 3) && (twoth < mythenMax)

Q = 2*pi*Y; % Scattering vector

% Not the Q=0 forward scattered reflection or Q > Qmax

if (h^2 + k^2 + l^2 ~= 0) && (hkl <= hklmax) && (1/Y > lambda/2)

% Determine multiplicity of each reflection depending

% on symmetry of the (hkl) values

if (k==0) && (l==0) % (h00) reflection, multiplicity = 6

M = 6;

end

if (k==h) && (l==0) % (hh0) reflection, multiplicity = 12

M = 12;

end

if (k==h) && (l==h) % (hhh) reflection, multiplicity = 8

M = 8;

end

if (k~=h) && (k~=0) && (l==0) % (hk0) reflection, multiplicity = 24

M = 24;

end

if (k~=h) && (k~=0) && (l==k) % (hkk) reflection, multiplicity = 24

M = 24;

end

if (k~=h) && (l~=h) && (l~=0) && (l~=k) && (k~=0) % (hkl) reflection, multiplicity = 48

M = 48;

end

% Now calculate form factors for silicon according

% to standard of four Gaussians plus a constant

fSi = a1*exp(-b1*(Q/(4*pi))^2) + ...

a2*exp(-b2*(Q/(4*pi))^2) + ...

a3*exp(-b3*(Q/(4*pi))^2) + ...

a4*exp(-b4*(Q/(4*pi))^2) + c;

FSi = fSi*(1+(-1)^(h+k) + (-1)^(k+l) + (-1)^(h+l))...

* (1 + (-1i)^(h+k+l)); % Structure factor for present h,k,l

if (FSi ~= 0)

sinth = Q*lambda/(4*pi); % sin(theta)

sin2th = sin(2*asin(sinth)); % sin(2theta)

LF = 1/(sinth*sin2th); % Lorentz factor

ISi = M*LF*FSi*conj(FSi); % Intensity of Bragg peak

twoth = 2*asind(sinth); % 2theta in degrees

nBPs = nBPs+1; % Increase number of registered Bragg peaks by 1

twothVector(nBPs) = twoth; % Add to vector array of 2theta positions of Bragg peaks

powderInt(nBPs) = ISi; % Add to vector array of intensities of Bragg peaks

% Add Bragg peak to powder pattern

powderPatty = powderPatty + ...

ISi*exp(-(powderPattx-twoth).^2/(2*sigma^2));

end % Determine intensities and 2theta values

end % Not direct beam and inaccessible peaks for which sin(theta) > 1

end % Only allowed peaks for silicon

end % l-loop

end % k-loop

end % h-loop

% -------------------------------------------------------------------

% Calculate x- and z-positions of curved plotted out diffraction pattern on

% the front face of the microstrip detector (@ y = -0.801)

RdataZero = 1.8*coneRad; % Radius of baseline with zero intensity

RdataMax = 1.2*coneRad; % Radius of maximum of diffraction data

Rdiff = RdataZero - RdataMax; % Span of data

powderPatty = powderPatty/max(powderPatty); % Normalize diffraction pattern to unity

Rdata = RdataZero - powderPatty*Rdiff; % Convert diffraction intensities to radii from (000)

Rdatax = -Rdata.*cosd(powderPattx); % x-values of curved data

Rdatay = 0*Rdatax-0.801; % y-values of curved data

Rdataz = -Rdata.*sind(powderPattx); % z-values of curved data

% -------------------------------------------------------------------

% Calculate x- and z-positions of standard Cartesian plot of diffraction

% pattern

xExtent = 20; % Width of plot in matlab units

xZero = -17; % x-origin in matlab units

zExtent = 10; % Height of plot in matlab units

zZero = -15.65; % z-origin in matlab units

% x, y, and z values for "normal" Cartesian plot of the powder pattern

powderCartx = powderPattx*xExtent/(max(powderPattx) - min(powderPattx))+xZero;

powderCarty = 0*powderPattx-0.801;

powderCartz = powderPatty*zExtent/(max(powderPatty) - min(powderPatty))+zZero;

% -----------------------------------------------------

% Determine surfs for Debye-Scherrer cones and strips

for kk = 1:nBPs % For all registered Bragg peaks

% Render entire cone

t = 0:1;

r = coneRad*t*sind(twothVector(kk));

[x,y,z] = cylinder(r,500);

surfArrayx(:,:,kk) = [x];

surfArrayy(:,:,kk) = [y];

surfArrayz(:,:,kk) = [z];

[x2,y2,z2] = cylinder(r,500);

for iii = 1:2

for jjj = 1:501

if (abs(y2(2,jjj))>0.4 || x2(2,jjj)>0)

x2(iii,jjj) = NaN;

y2(iii,jjj) = NaN;

z2(iii,jjj) = NaN;

end

stripArrayx(:,:,kk) = [x2];

stripArrayy(:,:,kk) = [y2];

stripArrayz(:,:,kk) = [z2];

end

% -------------------------------------------------------------------

% Set up movie parameters and graphics settings

vid = VideoWriter('PDexperimental.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(0, 'DefaultFigureRenderer', 'opengl');

set(gca,'linewidth',7);

myBlue = [0.4 0.44 0.73];

myCopper = [0.73 0.45 0.20];

myGold = [0.7969 0.6406 0.2383];

myBrass = [225 193 110]/255;

midGrey = [0.65 0.65 0.65];

lightGrey = [0.8 0.8 0.8];

lightGrey2 = [0.85 0.85 0.85];

darkGrey = [0.5 0.5 0.5];

lightBlue = [0.2 0.6 0.8];

lp = [1 -1 2];

lp2 = [1 0.1 0.5];

pos1 = [-0.1 -0.2 1.0 1.5];

maxAlpha = 0.5; % Maximum allowed opacity of Debye-Scherrer cone

maxPowderInt = max(powderInt); % Associated maximum signal intensity

moveSampleHolder = 3.3034; % Distance to move holder so middle of capillary is at (0,0,0)

angSum = 0;

frmNo = 0;

varAlp = 0.001;

% --------------

% Generate movie

subplot('position',pos1);

for ii = totAngleSteps % Loop according to vector containing change in angle for each new frame

frmNo = frmNo + 1; % Increase frame number by 1

angSum = angSum + ii; % Increase total rotational angle by value of ii

angNow = mod(angSum,360); % Just use values between 0 and 360 degrees

hold off;

newplot;

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

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

az = 40;

pol = 40;

set(gca,'view',[az,pol]);

axis equal

axis off

axis tight

hold on

% Dummy line to keep viewing field constant

plot3([-17 0],[-0.802 -0.802],[-20.1 -20.1],...

'color','w', 'LineWidth', 0.1);

% ------------------------------------------------------------

% Rotating capillary holder

[X,Y,Z] = cylinder(1,360);

baseRing1 = surf(X*0.95,Y*0.95,0.6*Z-moveSampleHolder,'FaceColor',...

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

rotate(baseRing1,[1 0 0],90,[0,0,0]);

rotate(baseRing1,[0 1 0],angNow,[0,0,0]);

baseRing2 = surf(X,Y,0.5*Z+0.05-moveSampleHolder,'FaceColor',...

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

rotate(baseRing2,[1 0 0],90,[0,0,0]);

rotate(baseRing2,[0 1 0],angNow,[0,0,0]);

theta = (0:360);

Xadd = 0.02*cosd(64*theta);

Yadd = 0.02*sind(64*theta);

X = X+Xadd; % Make ring nobbly

Y = Y+Yadd; % Make ring nobbly

fingerRing = surf(X,Y,0.4*Z+0.1-moveSampleHolder,'FaceColor',...

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

rotate(fingerRing,[1 0 0],90,[0,0,0]);

rotate(fingerRing,[0 1 0],angNow,[0,0,0]);

% Top circular surface of base ring

x = sind(theta);

y = cosd(theta);

z = 0.0*y;

topSurf1 = fill3(x,y,z+0.55-moveSampleHolder,lightGrey,'LineStyle','none',...

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

rotate(topSurf1,[1 0 0],90,[0,0,0]);

rotate(topSurf1,[0 1 0],angNow,[0,0,0]);

topSurf2 = fill3(0.95*x,0.95*y,z+0.6-moveSampleHolder,lightGrey,'LineStyle','none',...

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

rotate(topSurf2,[1 0 0],90,[0,0,0]);

rotate(topSurf2,[0 1 0],angNow,[0,0,0]);

% Lower translation stage

[V,F] = platonic_solid(2,3^0.5);

V(:,1) = 0.6*V(:,1);

V(:,2) = 0.6*V(:,2);

V(:,3) = 0.3*V(:,3)-moveSampleHolder+0.9;

stage1 = patch('Faces',F,'Vertices',V,'FaceColor',lightGrey,...

'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.5);

rotate(stage1,[1 0 0],90,[0,0,0]);

rotate(stage1,[0 1 0],angNow,[0,0,0]);

dovetail1 = plot3([-0.6 -0.3 -0.4 0.4 0.3 0.6,0.6 0.3 0.4 -0.4 -0.3 -0.6 -0.6],...

[0.6001 0.6001 0.6001 0.6001 0.6001 0.6001 -0.6001 -0.6001 -0.6001 -0.6001 -0.6001 -0.6001 0.6001],...

[0.8 0.8 0.9 0.9 0.8 0.8 0.8 0.8 0.9 0.9 0.8 0.8 0.8]+0.05-moveSampleHolder,'color',darkGrey, 'LineWidth',1.25);

rotate(dovetail1,[1 0 0],90,[0,0,0]);

rotate(dovetail1,[0 1 0],angNow,[0,0,0]);

[V,F] = platonic_solid(2,3^0.5);

V(:,1) = 0.05*V(:,1)-0.05;

V(:,2) = 0.1*V(:,2)-0.65;

V(:,3) = 0.05*V(:,3)-moveSampleHolder+0.8;

adjust1 = patch('Faces',F,'Vertices',V,'FaceColor',myBrass,...

'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.5);

rotate(adjust1,[1 0 0],90,[0,0,0]);

rotate(adjust1,[0 1 0],angNow,[0,0,0]);

% Upper translation stage

[V,F] = platonic_solid(2,3^0.5);

V(:,1) = 0.5*V(:,1);

V(:,2) = 0.5*V(:,2);

V(:,3) = 0.3*V(:,3)-moveSampleHolder+1.5;

stage2 = patch('Faces',F,'Vertices',V,'FaceColor',lightGrey2,...

'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.5);

rotate(stage2,[1 0 0],90,[0,0,0]);

rotate(stage2,[0 1 0],angNow,[0,0,0]);

dovetail2 = plot3([0.5001 0.5001 0.5001 0.5001 0.5001 0.5001 -0.5001 -0.5001 -0.5001 -0.5001 -0.5001 -0.5001 0.5001],...

[-0.5 -0.3 -0.4 0.4 0.3 0.5,0.5 0.3 0.4 -0.4 -0.3 -0.5 -0.5],...

[0.8 0.8 0.9 0.9 0.8 0.8 0.8 0.8 0.9 0.9 0.8 0.8 0.8]+0.65-moveSampleHolder,'color',darkGrey, 'LineWidth',1.25);

rotate(dovetail2,[1 0 0],90,[0,0,0]);

rotate(dovetail2,[0 1 0],angNow,[0,0,0]);

[V,F] = platonic_solid(2,3^0.5);

V(:,1) = 0.1*V(:,1)-0.55;

V(:,2) = 0.05*V(:,2)-0.025;

V(:,3) = 0.05*V(:,3)-moveSampleHolder+1.425;

adjust2 = patch('Faces',F,'Vertices',V,'FaceColor',myBrass,...

'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.5);

rotate(adjust2,[1 0 0],90,[0,0,0]);

rotate(adjust2,[0 1 0],angNow,[0,0,0]);

% Domed connector block

[X,Y,Z] = sphere(2000);

sphnot = 0.3;

% Only surviving part of sphere for X > xnot - a small slice of the sphere

X(abs(X)>sphnot) = NaN ; Y(abs(Y)>sphnot) = NaN ; Z(Z<sphnot) = NaN;

domeSurf = surf(X,Y,Z+1.25-moveSampleHolder,'FaceColor',lightGrey,'LineStyle','none',...

'FaceLighting','flat','DiffuseStrength',0.4);

rotate(domeSurf,[1 0 0],90,[0,0,0]);

rotate(domeSurf,[0 1 0],angNow,[0,0,0]);

radSide = (1^2 - sphnot^2)^0.5; % Radius of side wall of domed connector block

xx0 = -sphnot:0.01:sphnot; % Range of x-values for top curved part of side wall

angSide = asin(xx0); % Range of angles defining curved top side of side wall of domed connector

xx = [-sphnot xx0 sphnot -sphnot]; % Range of x-values for entire side wall

zSid = cos(angSide);

zSide = [0 zSid 0 0];

domeSide1 = fill3(xx,xx*0+sphnot,zSide+1.2034-moveSampleHolder,lightGrey,'LineStyle','none',...

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

rotate(domeSide1,[1 0 0],90,[0,0,0]);

rotate(domeSide1,[0 1 0],angNow,[0,0,0]);

domeSide2 = fill3(xx,xx*0+sphnot,zSide+1.2034-moveSampleHolder,lightGrey,'LineStyle','none',...

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

rotate(domeSide2,[0 0 1],180,[0,0,0]);

rotate(domeSide2,[1 0 0],90,[0,0,0]);

rotate(domeSide2,[0 1 0],angNow,[0,0,0]);

domeSide3 = fill3(xx,xx*0+sphnot,zSide+1.2034-moveSampleHolder,lightGrey,'LineStyle','none',...

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

rotate(domeSide3,[0 0 1],90,[0,0,0]);

rotate(domeSide3,[1 0 0],90,[0,0,0]);

rotate(domeSide3,[0 1 0],angNow,[0,0,0]);

domeSide4 = fill3(xx,xx*0+sphnot,zSide+1.2034-moveSampleHolder,lightGrey,'LineStyle','none',...

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

rotate(domeSide4,[0 0 1],270,[0,0,0]);

rotate(domeSide4,[1 0 0],90,[0,0,0]);

rotate(domeSide4,[0 1 0],angNow,[0,0,0]);

% Brass capillary holder ring

[X,Y,Z] = cylinder(1,360);

Xadd = 0.005*cosd(45*theta);

Yadd = 0.005*sind(45*theta);

X = X+Xadd;

Y = Y+Yadd;

fingerRing2 = surf(0.1*X,0.1*Y,0.1*Z+2.2034-moveSampleHolder,'FaceColor',...

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

rotate(fingerRing2,[1 0 0],90,[0,0,0]);

rotate(fingerRing2,[0 1 0],angNow,[0,0,0]);

% Top circular surface of base ring

x = 0.1*sind(theta);

y = 0.1*cosd(theta);

z = 0.0*y;

fingerRingTop = fill3(x,y,z+2.3024-moveSampleHolder,myBrass,'LineStyle','none',...

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

rotate(fingerRingTop,[1 0 0],90,[0,0,0]);

rotate(fingerRingTop,[0 1 0],angNow,[0,0,0]);

% Capillary glass - slightly off rotation axis (0.005) to lend a bit of

% wobbly verisimilitude to the animation

[X,Y,Z] = cylinder(1,360);

capillaryGlass = surf(0.035*X+0.005,0.035*Y,2*Z+2.3034-moveSampleHolder,...

lightGrey,'LineStyle','none','FaceLighting','gouraud',...

'FaceColor',lightGrey,'faceAlpha',0.2,'DiffuseStrength',1);

rotate(capillaryGlass,[1 0 0],90,[0,0,0]);

rotate(capillaryGlass,[0 1 0],angNow,[0,0,0]);

% Capillary powder - slightly off rotation axis (0.005) to lend a bit

% of wobbly verisimilitude to the animation

[X,Y,Z] = cylinder(1,360);

capillaryPowder = surf(0.014*X+0.005,0.014*Y,2*Z+2.3034-moveSampleHolder,...

'FaceColor',lightBlue,'LineStyle','none','FaceLighting','gouraud',...

'faceAlpha',1,'DiffuseStrength',1);

rotate(capillaryPowder,[1 0 0],90,[0,0,0]);

rotate(capillaryPowder,[0 1 0],angNow,[0,0,0]);

% -------------------------------------------------------------

% Create a cone of incident radiation

t = 0.35:0.05:1;

r = t/14;

[x,y,z] = cylinder(r,100);

if (frmNo <= sampleAccelFrames)

incRadAlp = 0.001; % Very weak but defines viewing region

elseif (frmNo > sampleAccelFrames) && (frmNo <= sampleAccelFrames+fadeInRadiation)

incRadAlp = (frmNo-sampleAccelFrames)/fadeInRadiation;

else

incRadAlp = 1;

end

sourceRad = surf(x/2,y,10*z,'FaceAlpha',incRadAlp,'FaceColor',...

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

rotate(sourceRad,[0 1 0],90,[0,0,0]);

% -------------------------------------------------------------

% Create set of Debye-Scherrer cones and strips of these that impinge

% on the microstrip detector

% Define Debye-Scherrer cone alpha (opacity) for all frames

fadeInConesStart = sampleAccelFrames;

fadeInConesFinish = fadeInConesStart+fadeInRadiation;

initialSitchStart = fadeInConesFinish;

initialSitchFinish = fadeInConesFinish+initialSitch;

fadeOutConesStart = initialSitchFinish;

fadeOutConesFinish = fadeOutConesStart+fadeOutFullFadeInDet;

if (frmNo <= sampleAccelFrames)

DSalp = 0.001; % Very weak but defines viewing region

elseif (frmNo > fadeInConesStart) && (frmNo <= fadeInConesFinish)

DSalp = (frmNo-fadeInConesStart)/fadeInRadiation; % Fade in DS cones

elseif (frmNo > initialSitchStart) && (frmNo <= initialSitchFinish)

DSalp = 1; % Full intensity cones

elseif (frmNo > fadeOutConesStart) && (frmNo <= fadeOutConesFinish)

DSalp = 1 - (frmNo-fadeOutConesStart)/fadeOutFullFadeInDet; % Fade out DS cones

else

DSalp = 0;

end

% Define Debye-Scherrer strip alpha (opacity) for all frames

if (frmNo <= fadeOutConesStart)

stripAlp = 0.001; % Very weak but defines viewing region

elseif (frmNo > fadeOutConesStart) && (frmNo <= fadeOutConesFinish)

stripAlp = (frmNo-fadeOutConesStart)/fadeOutFullFadeInDet; % Fade in DS strips

else

stripAlp = 1;

end

for kk = 1:nBPs % For all registered Bragg peaks

% Render entire cone

coneAlpha = maxAlpha*powderInt(kk)*DSalp/maxPowderInt;

sourceRad = surf(surfArrayx(:,:,kk),-surfArrayy(:,:,kk),...

coneRad*cosd(twothVector(kk))*surfArrayz(:,:,kk),'FaceAlpha',...

coneAlpha,'FaceColor',myGold,'LineStyle','none',...

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

rotate(sourceRad,[0 1 0],-90,[0,0,0]);

% Draw intersections of cones with microstrip detector brighter yellow

if (frmNo > sampleAccelFrames)

lineAlpha = incRadAlp*(powderInt(kk)/maxPowderInt)^0.7; % power of 0.7 to increase intensity of weak peaks

lineSig = plot3(stripArrayx(2,:,kk),-stripArrayy(2,:,kk),...

0.999*coneRad*cosd(twothVector(kk))*stripArrayz(2,:,kk),...

'color',[1 1 0 lineAlpha], 'LineWidth', 2);

rotate(lineSig,[0 1 0],-90,[0,0,0]);

end

% Render only that part of the cone falling on the detector strip

if (frmNo > fadeOutConesStart)

partConeAlpha = maxAlpha*powderInt(kk)*stripAlp/maxPowderInt;

stripRad = surf(stripArrayx(:,:,kk),-stripArrayy(:,:,kk),...

coneRad*cosd(twothVector(kk))*stripArrayz(:,:,kk)*1.001,...

'FaceAlpha',partConeAlpha,'FaceColor',myGold,'LineStyle',...

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

rotate(stripRad,[0 1 0],-90,[0,0,0]);

end

% -------------------------------------------------------------

% Create 1D microstrip detector

theta2 = (mythenMax-122:mythenMax+2);

x2 = cosd(theta2);

z2 = sind(theta2);

outCurvex2 = 2*coneRad*x2;

outCurvez2 = 2*coneRad*z2;

inCurvex2 = coneRad*x2;

inCurvez2 = coneRad*z2;

plateX = -[outCurvex2 flip(inCurvex2)];

plateZ = -[outCurvez2 flip(inCurvez2)];

plate1 = fill3(plateX,0*plateX+0.8,plateZ,lightGrey,'LineStyle','none',...

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

plate2 = fill3(plateX,0*plateX-0.8,plateZ,lightGrey,'LineStyle','none',...

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

curvePlateX = -[max(outCurvex2) inCurvex2 min(outCurvex2) min(outCurvex2)...

flip(inCurvex2) max(outCurvex2)];

curvePlateY = -[0.8 0*inCurvex2+0.8 0.8 -0.8 0*inCurvex2-0.8 -0.8];

curvePlateZ = -[outCurvez2(1) inCurvez2 outCurvez2(125) outCurvez2(125)...

flip(inCurvez2) outCurvez2(1)];

curvePlate = fill3(curvePlateX,curvePlateY,curvePlateZ,midGrey,'LineStyle',...

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

theta3 = (mythenMax-120:mythenMax);

x3 = cosd(theta3);

z3 = sind(theta3);

detSurfx3 = 0.999*coneRad*x3;

detSurfy3 = 0*detSurfx3;

detSurfz3 = 0.999*coneRad*z3;

detSurfX = -[detSurfx3 flip(detSurfx3)];

detSurfY = -[detSurfy3-0.4 flip(detSurfy3)+0.4];

detSurfZ = -[detSurfz3 flip(detSurfz3)];

detSurf = fill3(detSurfX,detSurfY,detSurfZ,'k','LineStyle','none',...

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

% --------------------------------------------------------------------

% Plot curved diffraction data on front surface of microstrip detector

% and then straighten it into a Cartesian plot

fadeInCurvedDataStart = fadeOutConesFinish;

fadeInCurvedDataFinish = fadeInCurvedDataStart+fadeInRadial;

if (frmNo <= fadeInCurvedDataStart)

curveAlp = 0;

elseif (frmNo > fadeInCurvedDataStart) && (frmNo <= fadeInCurvedDataFinish)

curveAlp = (frmNo - fadeInCurvedDataStart)/fadeInRadial;

else

curveAlp = 1;

end

% Straightening plot to Cartesian

radToLinStart = fadeInCurvedDataFinish+keepRadial;

radToLinFinish = radToLinStart+radToLin;

if (frmNo <= fadeInCurvedDataFinish+keepRadial)

plot3(Rdatax,Rdatay,Rdataz,...

'color',[0 0 0 curveAlp], 'LineWidth', 2.8);

elseif (frmNo > radToLinStart) && (frmNo <= radToLinFinish)

xCurveToLin = Rdatax - (Rdatax - powderCartx)*(frmNo-radToLinStart)/radToLin;

yCurveToLin = 0*Rdatax -0.801;

zCurveToLin = Rdataz - (Rdataz - powderCartz)*(frmNo-radToLinStart)/radToLin;

azvar = az - az*(frmNo - radToLinStart)/radToLin;

polvar = pol - pol*(frmNo - radToLinStart)/radToLin;

set(gca,'view',[azvar,polvar]);

plot3(xCurveToLin,yCurveToLin,zCurveToLin,...

'color','k', 'LineWidth', 2.8);

else % Last 90 frames showing Cartesian plot

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

plot3(powderCartx,powderCarty,powderCartz,...

'color','k', 'LineWidth', 2.8);

end

% -------------------------------------------------------------

% Store the frame

axis tight

axis off

frame = getframe(gcf);

writeVideo(vid,frame);

hold off

end

% Output the movie as an mpg file

close(vid);