% Program to plot out successive randomly oriented diffraction patterns
% recorded in an SSX experiment on lysozyme using a data file exptable.dat
% that was parsed from the original data file all.dat using the matlab
% routine parseDataFile.m

clear; close all;

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

% Open data file. This has the format:
%
% Chunk 1   h    k    l          I
% -38    7    1     -27.51
% -37    0    4      24.75  ...
% ...
% Chunk 2   h    k    l          I
% -32  -20   -4      47.25
% -30  -24    0      46.00  ...
% etc etc
%
% The program finds lines including the character 'h' and uses the data in
% between for each diffraction pattern

fid = fopen('exptable.dat','r');
S = textscan(fid,'%s','delimiter','\n');
SS = S{1};
fclose(fid);
% Get lines with 'h' in them - these are headers of data chunks
idx = find(contains(SS,'h')); % A list of line numbers for the headers
lengthFile = size(SS); % Number of lines in data file
ff = lengthFile(end,1);
numChunks = str2num(cell2mat(SS(ff)));

boxSize = 39.9999;
colormap jet;
% Quasitetragonal u.c. lattice constants of lysozyme in nm
abax = 7.86; % a = b
cax = 3.79; % c
cdab = cax/abax; % a/c

for i = 1:1000 % Number of diffraction patterns to plot. Gets REALLY slow for large i!!
    clear Pmat
    clear dumArray
    startInt = idx(i);
    endInt = idx(i+1);
    
    Pmat = zeros(endInt-startInt-1,4); % e.g. 5 to 500, 495 elements
    for j = endInt-startInt-1:-1:1
        Pmat(j,:) = str2num(cell2mat(SS(j+startInt))); % After much todo, at last a line of numbers
        % Determine if data point lies outside Q = 38 in a* and b* units
        h = Pmat(j,1);
        k = Pmat(j,2);
        l = Pmat(j,3);
        hkl = (h^2 + k^2 + (l/cdab)^2)^0.5;
        if (hkl >= 38) % Remove data above maximum plotted Q-value
            Pmat(j,:) = [];
        elseif (Pmat(j,4) <= 0) % Negative intensity data removed
            Pmat(j,:) = [];
        end
    end
    
    dumArray = [boxSize-2 0 0 0.00001; -(boxSize-2) 0 0 0.00001; ...
        0 boxSize-2 0 0.00001; 0 -(boxSize-2) 0 0.00001; ...
        0 0 (boxSize-2)*cdab 0.00001; 0 0 -(boxSize-2)*cdab 0.00001];
    circSize = 0.16*(0.000001+(Pmat(:,4)).^0.5);
    colorGrad = (Pmat(:,4)).^0.1; % Map colormap to intensity range of data
    if (i==1) % Start with an empty volume in reciprocal space
        scatter3(dumArray(:,1), dumArray(:,2), dumArray(:,3), dumArray(:,4), ...
            [0 0 0], 'MarkerFaceColor','Flat');
        xlim([-boxSize boxSize]);
        ylim([-boxSize boxSize]);
        zlim(cdab*[-boxSize boxSize]);
        % Drawn axes through origin. Matlab doesn't support this natively.
        % Comment out next seven lines if you want box axes
        line([0,0], [0,0], cdab*[-1.05*boxSize,1.05*boxSize], 'LineWidth', 2, 'Color', 'k');
        line([0,0], [-1.05*boxSize,1.05*boxSize], [0,0], 'LineWidth', 2, 'Color', 'k');
        line([-1.05*boxSize,1.05*boxSize], [0,0], [0,0], 'LineWidth', 2, 'Color', 'k');
        htext = text(-1.07*boxSize,0,0,'h','FontSize',14);
        ktext = text(0,-1.07*boxSize,0,'k','FontSize',14);
        ltext = text(-0.35,0,1.1*boxSize*cdab,'l','FontSize',14);
        
        axis square
        axis off % Comment out if you want box axes
        set(gca,'XTickLabel',[]);
        set(gca,'YTickLabel',[]);
        set(gca,'ZTickLabel',[]);
        set(gca, 'Projection','perspective');
        hold on
        
        pauseMov(vid,100);
    end
    scatter3(Pmat(:,1), Pmat(:,2), Pmat(:,3), circSize, colorGrad, 'MarkerFaceColor','Flat');
    LineHandle.NodeChildren.LineWidth = 0; % Set linewidth of outline of circles in scatter plot to zero
    
    cb = colorbar;
    set(cb,'position',[.35 .25 .01 .2],'TickLength',0.0); % Size and position of colorbar
    cb.Color = [1 1 1];
    
    % Comment next seven lines out for box axes
    line([0,0], [0,0], cdab*[-1.05*boxSize,1.05*boxSize], 'LineWidth', 2, 'Color', 'k');
    line([0,0], [-1.05*boxSize,1.05*boxSize], [0,0], 'LineWidth', 2, 'Color', 'k');
    line([-1.05*boxSize,1.05*boxSize], [0,0], [0,0], 'LineWidth', 2, 'Color', 'k');
    delete(htext); delete(ktext); delete(ltext)
    htext = text(-1.07*boxSize,0,0,'h','FontSize',14);
    ktext = text(0,-1.07*boxSize,0,'k','FontSize',14);
    ltext = text(-0.35,0,1.1*boxSize*cdab,'l','FontSize',14);
    
    axis square % Keep scaling 1:1:1 for x,y,z
    axis off % Comment out for box axes
    set(gca,'XTickLabel',[]);
    set(gca,'YTickLabel',[]);
    set(gca,'ZTickLabel',[]);
    
    if (i > 1)
        delete(hText);
    end
    str1 = 'Capture ';
    str2 = num2str(i,'% 3i');
    strTot = [str1,str2];
    hText = text(-1*boxSize, 0*boxSize, -0.82*boxSize/2,strTot,'FontSize',14);
    
    set(gca, 'Projection','perspective');
    hold on
    xlim([-boxSize boxSize]);
    ylim([-boxSize boxSize]);
    zlim(cdab*[-boxSize boxSize]);
    
    %Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    if (i <= 5)
        pauseMov(vid,15);
    end
    if (i >= 6) && (i <= 15)
        pauseMov(vid,7);
    end
    if (i >= 16) && (i <= 30)
        pauseMov(vid,4);
    end
    if (i >= 31) && (i <= 60)
        pauseMov(vid,3);
    end
    clear Pmat
    clear dumArray
end
delete(hText)
delete(htext)
delete(ktext)
delete(ltext)
[caz,cel] = view;
for i = caz:-caz/120:0 %120 default
    azimuth = i;
    elevation = cel*(i/caz);
    set(gca,'View',[azimuth,elevation]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

pauseMov(vid,45);
% Output the movie as an mpg file
close(vid);
