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 % 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

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);