Simulation of speckle patterns
Matlab code
% Program to generate movie of simulated XPCS speckle patterns for three
% different ensemble behaviours: (1) a single diffusing particle in an
% otherwise stationary; (2) an ensemble with increasing diffusion as
% temperature is increased linearly; (3) an ensemble beginning as a
% regular array crystal and melting; and (4) an ensemble of tumbling
% ellipsoids.
clear; close all
disp('Choose one of the following scenarios: Single diffusing particle (S),');
disp('Tumbling ellipsoids (E), Melting crystal (M),');
prompt = 'Diffusing sample with linearly increasing temperature (T): ';
str = input(prompt,'s');
if (str == 'S') || (str == 's')
vid = VideoWriter('singleParticleXPCS.mp4','MPEG-4');
frameMax = 300;
elseif (str == 'E') || (str == 'e')
vid = VideoWriter('tumblingEllipsoidsXPCS.mp4','MPEG-4');
frameMax = 300;
prompt = 'Ratio of long axis to short axes of ellipsoids: ';
bOverA = input(prompt);
elseif (str == 'M') || (str == 'm')
vid = VideoWriter('meltingCrystalXPCS.mp4','MPEG-4');
frameMax = 600;
elseif (str == 'T') || (str == 't')
vid = VideoWriter('diffusingEnsembleVtempXPCS.mp4','MPEG-4');
frameMax = 600;
end
pos3 = [0.68 0 0.32 1]; % Position of subplot of difference FT of ensemble between successive frames
prompt = 'Include Difference FT? (y/n)? ';
stryn = input(prompt,'s');
if (stryn == 'Y') || (stryn == 'y')
pos1 = [0 0 0.32 1]; % Position of subplot of ensemble of particles
pos2 = [0.34 0 0.32 1]; % Position of subplot of FT of ensemble
elseif (stryn == 'N') || (stryn == 'n')
pos1 = [0 0 0.5 1]; % Position of subplot of ensemble of particles
pos2 = [0.5 0 0.5 1]; % Position of subplot of FT of ensemble
end
vid.Quality = 100;
vid.FrameRate = 30; % 30 frame/s
open(vid);
figure('units','pixels','position',[0 0 4*4096 4*4096],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
% Ellipsoids require smaller relative minor axes or else it looks too
% crowded
if (str == 'E') || (str == 'e')
ra = 0.0125; % Radius of particles
R = ra*8; % Initial size of radius of ensemble
else
ra = 0.025; % Radius of particles
R = ra*4; % Initial size of radius of ensemble
end
xmax = R*2^4.31;
coords = zeros(2,1024); % Initialize radii and angles of particles' origins
% Generate set of coordinates for the centers of randomly distributed but
% more or less evenly distributed particles, or in the case of a melting
% crystal, a regular square array of circles
if (str ~= 'M') && (str ~= 'm') % Not option of melting crystal
for jj = 1:8
% Increase area by two to accommodate doubling of number of particles
R = R*sqrt(2);
if (jj == 1) || (jj == 2)
startval = jj;
else
startval = 2^(jj-2)+1;
end
endval = 2^(jj-1);
for ii = startval:endval % Double number of particles
% Place particles within annular ring with area equal to that of
% the circle within it
if (jj == 1)
theta = 0;
delr = 0;
coords(1,ii) = 0;
coords(2,ii) = 0;
else
theta = 2*pi*rand;
delr = R*(1 - 1/sqrt(2))*rand;
coords(1,ii) = (R/sqrt(2)+delr)*cos(theta);
coords(2,ii) = (R/sqrt(2)+delr)*sin(theta);
end
end
end
else
if (str == 'M') || (str == 'm') % Melting crystal
nUC = 0;
R = ra*4;
for jj = -8:8
for kk = -8:8
hold on
% Square lattice of identical circles
if (jj^2 + kk^2 <= 64)
nUC = nUC+1;
coords(1,nUC) = 16*R*jj/8;
coords(2,nUC) = 16*R*kk/8;
end % Inside radius
end % kk loop
end % jj loop
end % if M or m
endval = nUC; % Set frameMax to nUC
end
%__________________________________________________________________________
% Plot out circles or ellipsoids and record each frame
% Initial values for orientation of ellipsoids
phi1 = 0.5*pi*rand(1,endval); % Random phase between 0:90 degrees for out-of-plane rotation of ellipsoids
phi2 = pi*rand(1,endval); % Random phase between 0:180 degrees for in-plane rotation of ellipsoids
str2 = ' T_0';
th = 0:2*pi/50:2*pi; % Plotting angles of circles or ellipses
for frame = 1:frameMax
% Subplot of ensemble of particles
if (str == 'S') || (str == 's') % Single vibrating particle
pic1 = subplot('position',pos1);
newplot
hold off
difAmp = 2*ra*randn;
difAng = 2*pi*rand;
x_dif = difAmp*cos(difAng);
y_dif = difAmp*sin(difAng);
coords(1,1) = coords(1,1) + x_dif;
coords(2,1) = coords(2,1) + y_dif;
rdif = sqrt(coords(1,1)^2 + coords(2,1)^2); % Distance from centre
if (rdif > R) % Keep particle within circular boundary
coords(1,1) = coords(1,1) - x_dif;
coords(2,1) = coords(2,1) - y_dif;
end
% Plot out circles
for kk = 1:endval
hold on
x_circle = ra*cos(th) + coords(1,kk);
y_circle = ra*sin(th) + coords(2,kk);
plot(x_circle, y_circle);
if (kk == 1)
fill(x_circle, y_circle, 'r')
else
fill(x_circle, y_circle, 'k')
end
end
end
if (str == 'E') || (str == 'e') % Tumbling ellipsoids
pic1 = subplot('position',pos1);
newplot
hold off
x_ell = ra*sin(th);
y_ell = 3*ra*cos(th); % Ellipsoid with long axis three times the size of other two axes
% Plot out ellipses
for kk = 1:endval % endval = number of ellipses
hold on
difAmp = 0.5*ra*(1+randn);
difAng = 2*pi*rand;
x_dif = difAmp*cos(difAng);
y_dif = difAmp*sin(difAng);
coords(1,kk) = coords(1,kk) + x_dif;
coords(2,kk) = coords(2,kk) + y_dif;
rdif = sqrt(coords(1,kk)^2 + coords(2,kk)^2); % Distance from centre
if (rdif > R) % Keep particle within circular boundary
coords(1,kk) = coords(1,kk) - x_dif;
coords(2,kk) = coords(2,kk) - y_dif;
end
phi1(kk) = phi1(kk) + (pi/30)*randn; % New o-o-p angle
phi2(kk) = phi2(kk) + (pi/30)*randn; % New in-plane angle
y_ell = ra*cos(th)*(1+(bOverA-1)*(cos(phi1(kk)))^2); % Rotate ellipsoid out of plane
x_ell = ra*sin(th)*cos(phi2(kk)) + y_ell*sin(phi2(kk)); % and then in plane
y_ell = -ra*sin(th)*sin(phi2(kk)) + y_ell*cos(phi2(kk));
xvals = coords(1,kk)+x_ell;
yvals = coords(2,kk)+y_ell;
plot(xvals, yvals);
fill(xvals, yvals, 'k')
end
end
if (str == 'M') || (str == 'm') % Melting crystal
pic1 = subplot('position',pos1);
newplot
hold off
for kk = 1:endval % Draw each circle
hold on
% Square lattice begins to melt and diffuse into random
% positions
if (frame > 1)
difAmp = 0.125*ra*randn;
difAng = 2*pi*rand;
x_dif = difAmp*cos(difAng);
y_dif = difAmp*sin(difAng);
else % First frame
x_dif = 0;
y_dif = 0;
end
coords(1,kk) = coords(1,kk) + x_dif;
coords(2,kk) = coords(2,kk) + y_dif;
rdif = sqrt(coords(1,nUC)^2 + coords(2,nUC)^2); % Distance from centre
if (rdif > 16*R) % Keep particle within circular boundary
coords(1,kk) = coords(1,kk) - x_dif;
coords(2,kk) = coords(2,kk) - y_dif;
end
x_circle = ra*cos(th) + coords(1,kk);
y_circle = ra*sin(th) + coords(2,kk);
plot(x_circle, y_circle);
fill(x_circle, y_circle, 'k')
end
end
if (str == 'T') || (str == 't') % Diffusing ensemble with linearly increasing temperature
pic1 = subplot('position',pos1);
newplot
hold off
for kk = 1:endval
difAmp = 0.01*sqrt(frame)*ra*randn;
difAng = 2*pi*rand;
x_dif = difAmp*cos(difAng);
y_dif = difAmp*sin(difAng);
hold on
% Randomly positioned diffusing identical circles
coords(1,kk) = coords(1,kk) + x_dif;
coords(2,kk) = coords(2,kk) + y_dif;
rdif = sqrt(coords(1,kk)^2 + coords(2,kk)^2); % Distance from centre
if (rdif > R) % Keep particle within circular boundary
coords(1,kk) = coords(1,kk) - x_dif;
coords(2,kk) = coords(2,kk) - y_dif;
end
x_circle = ra*cos(th) + coords(1,kk);
y_circle = ra*sin(th) + coords(2,kk);
plot(x_circle, y_circle);
fill(x_circle, y_circle, 'k')
end
end
% Add white circles at corners to stop saved image shrinking to size of
% present ensemble
for aa = -1:2:1
for bb = -1:2:1
x_circle = ra * cos(th) + aa*xmax;
y_circle = ra * sin(th) + bb*xmax;
plot(x_circle, y_circle,'w*','Linewidth',0.0001);
pp = fill(x_circle, y_circle, 'w');
pp.LineWidth = 0.0001;
pp.EdgeColor = [0.99 0.99 0.99];
end
end
axis square
xlim([-xmax xmax])
ylim([-xmax xmax])
axis off
filename = 'circles.png';
myAxes=findobj(pic1,'Type','Axes');
exportgraphics(myAxes,filename,'Resolution',300);
str1 = num2str(frame,'%3i');
if (str == 'T') || (str == 't')
if (frame == 1)
str1 = ' ';
end
strTot = [str1,str2];
hText = text(1.25*R,-1.25*R,strTot, 'FontSize',25,'VerticalAlignment',...
'bottom','HorizontalAlignment','right');
end
hold off
%______________________________________________________________________
% Subplot of FT of ensemble
subplot('position',pos2);
newplot
I = imread(filename);
hold off
I1 = double(I);
I1=imcomplement(rgb2gray(I1));
F = fft2(I1);
F = abs(fftshift(F));
if (frame == 1)
Fprev = abs(F); % Previous FT
end
minInt = min(min((abs(F)).^0.32)); % Power 0.4 compresses intensity scale
maxInt = max(max(abs(F))); % Find central maximum
[xxx,yyy]=find(F==maxInt);
xsize = round(size(F,1)/10); % 10% of linear dimension of FT image
Fred = F(xxx-xsize:xxx+xsize,yyy-xsize:yyy+xsize); % Central 20% x 20% of FT
maxInt = max(max((abs(Fred)).^0.32)); % Power 0.4 compresses intensity scale
imshow(((abs(Fred)).^0.32),[minInt,maxInt]);
greenBlack = customcolormap([0 1], {'#00ff00','#000000'}); % Similar to laser green speckle
colormap(gca,greenBlack);
%______________________________________________________________________
if (stryn == 'Y') || (stryn == 'y') % Include difference FT
subplot('position',pos3); % Subplot of relative differences between successive FTs
newplot
DiffFT = 2*(F - Fprev)./(F + Fprev); % Varies between limits of +/-2
DiffFT = DiffFT(xxx-xsize:xxx+xsize,yyy-xsize:yyy+xsize); % Central 20% x 20% of FT
maxInt = max(max((abs(DiffFT))));
imshow((((DiffFT))),[-2,2]);
RWB = customcolormap([0 0.5 1], {'#0000ff','#ffffff','#ff0000'}); % Red-White-Blue
colormap(gca,RWB);
cb = colorbar;
cb.Position = [0.975,0.25,0.01,0.16];
caxis([-2 2]);
set(cb,'TickLabels',{'-2','0','2'},'Ticks',[-2,0,2],...
'FontSize',14);
set(cb,'YAxisLocation','left')
end
frame2 = getframe(gcf);
writeVideo(vid,frame2);
Fprev = F; % Set present FT to be previous FT in next loop cycle
hold off
end
close(vid);