Reciprocal-space 3D plot of SSX data for first 1000 patterns
Matlab code
% 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);
% Parsing program of huge Crystfel data file to extract only h k l
% and intensity up to a certain h k l, ignoring negative intensities.
% Note: this program will produced a parsed file with a size in the range
% of 100s of MB!! It takes a couple of minutes to generate on my 2018
% macbook pro, too!
clear; close all;
Str = fileread('all.dat');
CStr = strsplit(Str, newline);
N = size(CStr,2);
fileID = fopen('exptable.dat','w');
fileChunk=0;
printToFile = 0;
for i = 1:N
% Find start of data chunk
lineId = contains(CStr(i), 'h k l'); % Logical identifier of lines that start new data chunk
if (lineId == logical(1)) % i.e. YES! this line contains 'h k l'!!
fileChunk = fileChunk+1;
%CStr(i)
fprintf(fileID,'Chunk %i',fileChunk);
printToFile = 1; % Allows next line to be printed to file
end
% Find end of data chunk
lineId2 = contains(CStr(i), 'End of reflections'); % Line below the end of latest data chunk
if (lineId2 == logical(1)) % i.e. YES! this line contains 'End of reflections'!!
%CStr(i-1)
printToFile = 0; % Stops next line to be printed to file
end
% Writing to file part of loop
if (printToFile == 1)
thisLine = cell2mat(CStr(i));
fprintf(fileID,'%28.28s \n',thisLine);
end
end
fprintf(fileID,'%i \n',fileChunk);fileChunk % Number of diffraction patterns
fclose(fileID);