% 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; and (3) an ensemble beginning as a
% regular array crystal and melting

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