Simulation of speckle patterns 

Matlab code 

% Program to generate movie of scattering pattern of increasing number of

% particles but at a roughly constant particle density whereby 

% interference between scattering from individual particles generates 

% speckle, whose characteristic separation is inversely proportional to 

% the linear size of the envelope of the ensemble of the particles. Three

% options are possible: identical circles, variable-radius circles, or

% variable size polygons. 

 

clear; close all

 

disp('Choose one of the following ensembles: Identical circles (I),');

prompt = 'variable-radius circles (V), Random polygons (P): ';

str = input(prompt,'s');

if (str == 'I') || (str == 'i')

    vid = VideoWriter('identicalCircleSpeckle.mp4','MPEG-4');

elseif (str == 'V') || (str == 'v')

    vid = VideoWriter('variableCircleSpeckle.mp4','MPEG-4');

else

    vid = VideoWriter('polygonSpeckle.mp4','MPEG-4');

end

vid.Quality = 100;

vid.FrameRate = 1; % 1 frame/s 

open(vid);

figure('units','pixels','position',[0 0 4*4096 4*4096],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);

% Plotting areas of two components of video

pos1 = [0 0 0.5 1]; % Ensemble of particles 

pos2 = [0.5 0 0.5 1]; % FT of ensemble

ra = 0.04; % Radius of particles

R = ra*2.5; % Initial size of radius of ensemble 

xmax = R*2^4.52;

coords = zeros(2,1024); % Initialize radii and angles of circle origins

for jj = 1:9 

    % Subplot of ensemble of particles 

    pic1 = subplot('position',pos1); 

    % 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 

        if (jj == 1)

            theta = 0;

            delr = 0;

        else

            theta = 2*pi*rand;

            delr = R*(1 - 1/sqrt(2))*rand;

        end 

        % Place particles within annular ring with area equal to circle

        % within it 

        coords(1,ii) = (R/sqrt(2)+delr)*cos(theta); 

        coords(2,ii) = (R/sqrt(2)+delr)*sin(theta);

    end

    hold off

    for kk = 1:endval

        hold on

        if (str == 'P') || (str == 'p') % Polygons and circles 

            polyVal = 0.2 + rand*0.8; % Determines number of sides of polygon 

            polyVarSize = 0.5+rand; % Determines size of polygon 

            if (polyVal >= 0.3) % Not a circle

                sides = round(10*polyVal);

                th = 0:2*pi/sides:2*pi;

                randAng = 2*pi*rand; % Determines angular orientation of polygon 

                x_circle = ra*polyVarSize*cos(th+randAng) + coords(1,kk);

                y_circle = ra*polyVarSize*sin(th+randAng) + coords(2,kk);

                plot(x_circle, y_circle);

                fill(x_circle, y_circle, 'k')

            else % Circle 

                th = 0:2*pi/50:2*pi;

                x_circle = ra*polyVarSize*cos(th) + coords(1,kk);

                y_circle = ra*polyVarSize*sin(th) + coords(2,kk);

                plot(x_circle, y_circle);

                fill(x_circle, y_circle, 'k')

            end

        elseif (str == 'V') || (str == 'v') % Variable-radius circles 

            th = 0:2*pi/50:2*pi;

            VarSize = 0.5+rand; % Radius between 0.5 and 1.5 of average 

            x_circle = ra*VarSize*cos(th) + coords(1,kk);

            y_circle = ra*VarSize*sin(th) + coords(2,kk);

            plot(x_circle, y_circle);

            fill(x_circle, y_circle, 'k')

        else % Identical circles 

            th = 0:2*pi/50:2*pi;

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

    

    % Subplot of FT of ensemble 

    subplot('position',pos2); 

    I = imread(filename);

    hold off

    I1 = double(I);

    I1=imcomplement(rgb2gray(I1));

    F = fft2(I1);

    F = fftshift(F);

    minInt = min(min((abs(F)).^0.4)); % Power 0.4 compresses intensity scale

    maxInt = max(max((abs(F)).^0.4));

    xsize = round(0.1*size(F,1)); % 10% of linear dimension of FT image 

    F = F(4*xsize:6*xsize,4*xsize:6*xsize); % Central 20% x 20% of FT 

    imshow(((abs(F)).^0.4),[minInt,maxInt/(1+jj/5)]); % Scale varies with jj to keep nice intensity range

    greenBlack = customcolormap([0 1], {'#00ff00','#000000'}); % Similar to laser green speckle 

    colormap(greenBlack);

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

close(vid)