% 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.2; % degrees increase/frame during sample angular acceleration 
topSpeed = 15; % 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 

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 = -20; % 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; 

% -------------------------------------------------------------------
% 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(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;
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');
    set(gca,'view',[40,40]);
    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(45*theta);
    Yadd = 0.02*sind(45*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','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  
        t = 0:1;
        r = coneRad*t*sind(twothVector(kk));
        [x,y,z] = cylinder(r,500);
        coneAlpha = maxAlpha*powderInt(kk)*DSalp/maxPowderInt; 
        sourceRad = surf(x,-y,coneRad*cosd(twothVector(kk))*z,'FaceAlpha',...
            coneAlpha,'FaceColor',...
            myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(sourceRad,[0 1 0],-90,[0,0,0]);
        
        % Render only that part of the cone falling on the detector strip 
        t = 0:1;
        r = coneRad*t*sind(twothVector(kk));
        [x,y,z] = cylinder(r,500);
        for iii = 1:2
            for jjj = 1:501
                if (abs(y(2,jjj))>0.4 || x(2,jjj)>0)
                    x(iii,jjj) = NaN;
                    y(iii,jjj) = NaN;
                    z(iii,jjj) = NaN;
                end
            end
        end
        partConeAlpha = maxAlpha*powderInt(kk)*stripAlp/maxPowderInt; 
        sourceRad2 = surf(x,-y,coneRad*cosd(twothVector(kk))*z*1.001,'FaceAlpha',...
            partConeAlpha,'FaceColor',...
            myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(sourceRad2,[0 1 0],-90,[0,0,0]); 
        % Draw intersections of cones with microstrip detector brighter yellow
        lineAlpha = incRadAlp*(powderInt(kk)/maxPowderInt)^0.7; % power of 0.7 to increase intensity of weak peaks 
        lineSig = plot3(x(2,:),-y(2,:),0.999*coneRad*cosd(twothVector(kk))*z(2,:),...
            'color',[1 1 0 lineAlpha], 'LineWidth', 2);
        rotate(lineSig,[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 
    fadeInCurvedDataStart = fadeOutConesFinish;
    fadeInCurvedDataFinish = fadeInCurvedDataStart+fadeInRadial;
    if (frmNo <= fadeInCurvedDataStart) 
        curveAlp = 0; 
    elseif (frmNo > fadeInCurvedDataStart) && (frmNo <= fadeInCurvedDataFinish)
        curveAlp = (frmNo - fadeInCurvedDataStart)/fadeInRadial;
    else 
        curveAlp = 1; 
    end 
    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; 
        plot3(xCurveToLin,yCurveToLin,zCurveToLin,...
            'color','k', 'LineWidth', 2.8);
    else % Last 90 frames showing Cartesian plot 
        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);
