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