% Cartoon of an angle-scanned XPD experiment. Data from Matthias Muntwiler
% and Martin Heinrich, PEARL beamline, Swiss Light Source 
close all; clear

vid = VideoWriter('xpdExpt.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);
LL = 1.05;
[x,y,z] = sphere(70); % Sphere used for atoms
[X,Y,Z] = sphere(90); % Shaded sphere on which the XPD data is plotted
X = X(46:end,:);       % Keep top half
Y = Y(46:end,:);       % Keep top half
Z = Z(46:end,:);       % Keep top half
% Shrink size of sphere a bit to stop it popping through the data pixels now
% and again
X = 0.998*X;
Y = 0.998*Y;
Z = 0.998*Z;

GeTe = importdata('GeTe_Te4s_626eV.dat'); % XPD data
GeTe = flip(GeTe,1);

% Steps in polar direction are constant and equal to 1 degrees. Upper and
% lower limits of each data pixel given by theta1 and theta2
theta = GeTe(:,2); % Polar angles
theta1 = theta-0.55; % 0.55 instead of 0.5 to ensure overlap
theta2 = theta+0.55;
phi = 360-GeTe(:,3); % Azimuthal angles
len = size(phi,1); % Number of data points
% Convert phi so it increases only and doesn't reset to zero in next
% azimuthal scan. This is needed to calculated phi1 and phi2, the left and
% right limits of the data square
kk = 1;
for ii = 1:len-1
    if (phi(ii+1) < phi(ii)) % New azimuthal scan detected
        phi(ii+1:len) = phi(ii+1:len)+360; % Increase phi from hereon by 360 deg
        phi0deg(kk) = ii; % Vector of data-point indices for which new phi rotation begins 
        kk = kk+1;
    end
end
phinxt = circshift(phi,-1); % Next phi value brought back to same index as present phi value
phiprv = circshift(phi,+1); % Previous phi value brought forward to same index as present phi value
% Left and right limits of data square determined by phi1 and phi2
phi1 = 0.1+(phinxt+phi)/2;
phi1(len) = 2*phi1(len-1)-phi1(len-2);
phi2 = circshift(phi1,+1)-0.1; % What was, is now. What is now, will be next
phi2(1) = 2*phi2(2)-phi2(3); % Subtract 0.2 to remove small gap at start/end of phi scans

inty = GeTe(:,1); % Intensity values of data at theta, phi

% Convert from spherical polar coordinates to Cartesian coordinates
X1 = sind(theta1).*cosd(phi1);
X2 = sind(theta1).*cosd(phi2);
X3 = sind(theta2).*cosd(phi2);
X4 = sind(theta2).*cosd(phi1);
Y1 = sind(theta1).*sind(phi1);
Y2 = sind(theta1).*sind(phi2);
Y3 = sind(theta2).*sind(phi2);
Y4 = sind(theta2).*sind(phi1);
Z1 = cosd(theta1);
Z2 = cosd(theta1);
Z3 = cosd(theta2);
Z4 = cosd(theta2);

% Get structural data of crystalline sample unit cell. Contains x, y, z
% coordinates of each atom, its atomic number, and its radius (ionic or
% covalent, whichever is more appropriate)

sampleData = importdata('GeTeCoords.dat');
noAtoms = size(sampleData,1);
C = zeros(noAtoms,3);
GeCol = [102 143 143]/256;
TeCol = [212 122 0]/256;
for i = 1:noAtoms
    if (sampleData(i,1) <= 1.28)
        C(i,1:3) = GeCol;
    else
        C(i,1:3) = TeCol;
    end
end

hold off

% Color index for "rgb" to determine color of each data square based on its
% intensity
rgb = customcolormap([0 0.5 1],{'#ffff00','#ff0000','#000000'});
Call = zeros(len,3);
colIndex = round(256*(inty-min(inty))/(max(inty)-min(inty)));
for ii = 1:len
    if (colIndex(ii) == 0)
        colIndex(ii) = 1;
    end
    Call(ii,:) = rgb(colIndex(ii),:);
end
colormap(Call);

% OK. All prep stuff done. Now for the plotting...
pt1 = 2:2:phi0deg(3);
pt2 = phi0deg(4:end);
allphiIndex = [pt1 pt2 len];
% And we're off! On the way-out Wacky Races!
for jj = allphiIndex % Data point loop
    newplot
    col = 1:jj;
    
    tic % Remove, also toc, if you're not obsessed with how terribly
    % long it takes to render this video. Planning to buy a new mac
    % soon...
    for ii = 1:jj % Plot each pixel recorded so far
        xcoords = [X1(ii) X2(ii) X3(ii) X4(ii)];
        ycoords = [Y1(ii) Y2(ii) Y3(ii) Y4(ii)];
        zcoords = [Z1(ii) Z2(ii) Z3(ii) Z4(ii)];
        datPix = fill3(xcoords,ycoords,zcoords,Call(ii,:),...
            'linestyle','none','FaceAlpha',1);
        if (jj <= phi0deg(3)) % Plot out phi rotation for first three full rotations only 
            rotate(datPix,[0 0 1],-phi1(jj),[0 0 0]) % Rotate around z-axis by phi
        end
        rotate(datPix,[0 1 0],90-theta(jj),[0 0 0]) % Rotate around y-axis by theta
        hold on
    end
    
    % Plot atoms of crystalline sample
    scaleFac = 25;
    for i = 1:noAtoms
        rad = sampleData(i,1)/scaleFac;
        a = sampleData(i,2)/scaleFac;
        b = sampleData(i,3)/scaleFac;
        c = sampleData(i,4)/scaleFac;
        atomCol = C(i,:);
        % Create concatenated surface image of the unit cell
        sample = surf(rad*x+a,rad*y+b,rad*z+c,...
            'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
            'none','FaceLighting','flat','AmbientStrength',0.2,...
            'DiffuseStrength',0.5,'Visible','off');
        hold on
        if (jj <= phi0deg(3))
            rotate(sample,[0 0 1],-phi1(jj),[0 0 0]) % Rotate around z-axis by phi
        end
        rotate(sample,[0 1 0],90-theta(jj),[0 0 0]) % Rotate around y-axis by theta
        set(sample,'Visible','on')
    end
    % Axis to entrance of CHA 
    plot3([0 2],[0 0],[0 0],'LineWidth',1.7,'LineStyle','-.','Color','y') 
    % Hemisphere axis, normal to sample surface 
    hsAxis = plot3([0 0],[0 0],[0 2],'LineWidth',1.7,'LineStyle','-.','Color','k'); 
    hemisphere = surf(X,Y,Z,'FaceAlpha',0.25,'FaceColor','k','LineWidth',1,'EdgeColor',[0.7 0.7 0.7],'FaceLighting','flat',...
        'DiffuseStrength',1);
    material(hemisphere,'shiny');
    if (jj <= phi0deg(3))
        rotate(hemisphere,[0 0 1],-phi1(jj),[0 0 0]) % Rotate around z-axis by phi
    end
    rotate(hemisphere,[0 1 0],90-theta(jj),[0 0 0]) % Rotate hemisphere around y-axis by theta
    rotate(hsAxis,[0 1 0],90-theta(jj),[0 0 0]) % Rotate axis around y-axis by theta
    
    lp = [0 -1 5];
    lp2 = [7 2 2];
    light('Position',lp,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    set(gca, 'Projection','perspective');
    set(gca,'View',[116,0]);
    axis equal
    axis off
    
    xlim([-LL LL+0.2]);
    ylim([-LL LL]);
    zlim([-1 LL]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
    if (jj > phi0deg(3))
        pauseMov(vid,1) % Slow down filling in of data a bit
    end
    toc
    hold off
end % End data point loop
close(vid)