Order_4_carpet_anisotropic.m

% Program made by Gabriel María Ybarra Marcaida 
% under Luis Vega González's supervision.

% Memory clear (just in case)
clear
clearvars;
clearvars -global;

% ------ DEFINITION OF THE VARIABLES ------

% PHYSICAL PARAMETERS

A = 1;                      % Wave's amplitude
c = 1;                      % Speed of light
lambda = 0.1;               % Wavelength
omega = 2*pi*c / lambda;    % Angular frequency
w = 0.2*lambda;             % Grating width.
d = 1;                      % Grating separation.
zt = 2*d^2/lambda;          % Talbot distance.


% CANVAS PARAMETERS

X_max = d/2;
Z_max = zt;

extra_steps = 1.52;         % The total number of steps is zt*extra_steps/h_z

% DISCRETIZATION PARAMETERS

N_x = 20;                   % Space discretization, with h_x = λ/N_x.
N_z = 0.1*N_x^2;            % Space discretization, with h_z = λ/N_z.

M = sqrt(2)*N_z;                        % Time discretization, with τ = f/M.
r = ( 1/M / (1/N_x + 1/N_z) )^2;        % Must be smaller than 1 if we want a stable algorithm.
a_discrete = ceil( w/2/lambda*N_x );    % Grating's discrete width.

tau = lambda/c/M; % Temporal step, τ = f/M.
h_x = lambda/N_x; % Spatial X step (every matrix element represents an area of hxh). hX = λ/N.
h_z = lambda/N_z; % Spatial Z step (every matrix element represents an area of hxh). hZ = λ/N.

n_x = ceil( X_max/h_x ); % Number of points in the x direction.
n_z = ceil( Z_max/h_z ); % Number of points in the z direction.



% OTHER VARIABLES

t = 0;                                      % Current time
x = double((-3*X_max:double(h_x):3*X_max)); % Matrix with the x-axis positions
z = double((0:double(h_z):Z_max));          % Matrix with the z-axis positions


tstep = 1;  % tstep counter
modulo = 3; % Number of tsteps between photos






% -------- BEGINING OF THE MAIN PROGRAM -----------


% We check that r is smaller than 1/2
if r > 1
    disp("r = " + num2str(r) + " which is less than 1/2!");
    disp("The algorithm isn't stable!");
    return;
end



% Definition of the inicial matrices
u_ahora = zeros(n_x + 5,floor(n_z*extra_steps) + 4);
u_despues = u_ahora;
u2 = zeros(n_x+1,n_z+1);

% CALCULATION OF U AT N=1
t = tau;

% WE SET THE SOURCE
u_despues(3:a_discrete+3,1) = A*sin(omega*(t-tau));
u_despues(3:a_discrete+3,2) = A*sin(omega*t);

% WE IMPOSE PERIODICITY
u_despues(1,:) = u_despues(5,:);
u_despues(2,:) = u_despues(4,:);

u_despues(end-4,:) = u_despues(end,:);
u_despues(end-3,:) = u_despues(end-1,:);

% We update the matrices
u_antes = u_ahora;
u_ahora = u_despues;

for j=1:tstep
    u2(:,j) = u_ahora(3:end-2,j).^2;
end


% WE GENERATE THE IMAGE

fig = figure;
fig.WindowState = 'maximized';

set(gca,'FontSize',20)


colormap('turbo');  % set colormap
colorbar;           % set colorbar

foldername = ['images_d=',num2str(d),'_w=',num2str(w),'_l=',num2str(lambda),'_N=',num2str(N_z)];
destdirectory = fullfile(pwd,foldername);
mkdir(destdirectory);

generate_image(x,z,tstep,u2,A,fig,destdirectory);



% Calculamos el estado del sistema en el tiempo t = k*l*tau
    
    for tstep = 2:n_z
        disp("Step " + num2str(tstep) + " of " + num2str(n_z*extra_steps) + ".");
        t = t + tau;
        
        
        % WE SET THE SOURCE
        u_despues(3:a_discrete+3,1) = A*sin(omega*(t-tau));
        u_despues(3:a_discrete+3,2) = A*sin(omega*t);
        
        
        % We update every point
        for m = 3:tstep
            for l = 3:n_x+3
                u_despues(l,m) = 2*u_ahora(l,m) - u_antes(l,m) +...
                    (c*tau)^4 /12 * laplacian_squared(u_ahora,l,m,h_x,h_z) +...
                    (c*tau)^2 * laplacian(u_ahora,l,m,h_x,h_z);
            end
        end
        
        % WE IMPOSE PERIODICITY
        u_despues(1,1:tstep) = u_despues(5,1:tstep);
        u_despues(2,1:tstep) = u_despues(4,1:tstep);

        u_despues(end,1:tstep) = u_despues(end-4,1:tstep);
        u_despues(end-1,1:tstep) = u_despues(end-3,1:tstep);
        
        
        % We update the matrices
        u_antes = u_ahora;
        u_ahora = u_despues;
        

        
        % We generate the image
        if mod(tstep,modulo) == 1 || modulo==1

            % We compute the intensity
            for j=1:n_z+1
                u2(:,j) = u_ahora(3:end-2,j+1).^2;
            end

            generate_image(x,z,tstep,u2,A,fig,destdirectory);
        end
        
    end

    for tstep = n_z+1:floor(n_z*extra_steps)
        disp("Step " + num2str(tstep) + " of " + num2str(n_z*extra_steps) + ".");
        t = t + tau;
        
        
        % WE SET THE SOURCE
        u_despues(3:a_discrete+3,1) = A*sin(omega*(t-tau));
        u_despues(3:a_discrete+3,2) = A*sin(omega*t);
        
        
        % We update every point
        for m = 3:tstep
            for l = 3:n_x+3
                u_despues(l,m) = 2*u_ahora(l,m) - u_antes(l,m) +...
                    (c*tau)^4 /12 * laplacian_squared(u_ahora,l,m,h_x,h_z) +...
                    (c*tau)^2 * laplacian(u_ahora,l,m,h_x,h_z);
            end
        end
        
        % WE IMPOSE PERIODICITY
        u_despues(1,1:tstep) = u_despues(5,1:tstep);
        u_despues(2,1:tstep) = u_despues(4,1:tstep);

        u_despues(end,1:tstep) = u_despues(end-4,1:tstep);
        u_despues(end-1,1:tstep) = u_despues(end-3,1:tstep);
        
        
        % We update the matrices
        u_antes = u_ahora;
        u_ahora = u_despues;
        

        % We generate the image
        if mod(tstep,modulo) == 1 || modulo==1

            % We compute the intensity
            for j=1:n_z+1
                u2(:,j) = u_ahora(3:end-2,j+1).^2;
            end

            generate_image(x,z,tstep,u2,A,fig,destdirectory);
        end

    end

    

% We generate the video
generate_video(tstep,modulo,destdirectory);



fprintf("\n\nProgram made by Gabriel María Ybarra Marcaida under the supervision of Luis Vega González.\n");
return;




% --------- FUNCTIONS ---------


% Computes the laplacian of u at (l,m)
% A fourth order scheme is used
function lapl = laplacian(u,l,m,h_x,h_z)

    lapl = -5/2*u(l,m) / h_x^2 - 5/2*u(l,m) / h_z^2 + ...
           4*( u(l+1,m) + u(l-1,m) )/3 / h_x^2 + ...
           4*( u(l,m+1) + u(l,m-1) )/3 / h_z^2 + ...
           - ( u(l+2,m) + u(l-2,m) )/12 / h_x^2 + ...
           - ( u(l,m+2) + u(l,m-2) )/12 / h_z^2;
    
end

% Computes the squared laplacian of u at (l,m)
% A second order scheme is used
function lapl2 = laplacian_squared(u,l,m,h_x,h_z)
    
    lapl2 = ( u(l-2,m) - 4*u(l-1,m) + 6*u(l,m) - 4*u(l+1,m) + u(l+2,m) ) / h_x^4 + ...
            ( u(l,m-2) - 4*u(l,m-1) + 6*u(l,m) - 4*u(l,m+1) + u(l,m+2) ) / h_z^4 + ...
            2*( 4*u(l,m) + ( u(l+1,m+1) + u(l+1,m-1) + u(l-1,m+1) + u(l-1,m-1) ) - 2*( u(l+1,m) + u(l-1,m) + u(l,m+1) + u(l,m-1) ) ) / h_x^2 / h_z^2;

end


% Generates the image

function generate_image(x,z,tstep,u2,A,fig,destdirectory)

    imagen = [flip(u2(3:end-3,:)); u2(3:end-2,:); flip(u2(4:end-3,:)); u2(3:end-2,:); flip(u2(4:end-3,:)); u2(3:end-2,:)];

    imagesc(x,z,imagen, [0 A^2]);   % draw image and scale colormap to values range
    colorbar;                       % set colorbar

    set(gca,'FontSize',20)

    Z = length(z)-1;
    yticks([0 z(floor(Z/6)) z(floor(Z/3)) z(floor(Z/2)) z(floor(2*Z/3)) z(floor(5*Z/6)) z(Z)])
    yticklabels({'-3/2 d','- d','-1/2 d','0','1/2 d','d','3/2 d'})
    
    
    X = length(x)-1;
    xticks([x(1) x(floor(X/4)) x(floor(X/2)) x(floor(3*X/4)) x(X)])
    xticklabels({'0','1/4 Z_T','1/2 Z_T','3/4 Z_T','Z_T'})

    drawnow
    frame = getframe(fig);
    im = frame2im(frame);
    thisimage = ['img_',num2str(tstep),'.png'];
    fulldestination = fullfile(destdirectory,thisimage);
    imwrite(im,fulldestination);
end

function generate_video(tstep,modulo,destdirectory)
    % Create a video writer object
    writerObj = VideoWriter(fullfile(destdirectory,'aa'));

    % Set frame rate
    writerObj.FrameRate = 30;

    % Open video writer object and write frames sequentially
    open(writerObj)

    for i = 1:modulo:tstep                   % Some number of frames
         % Read frame
         frame = sprintf('img_%d.png', i);
         frame = fullfile(destdirectory, frame);
         input = imread(frame);

         % Write frame now
         writeVideo(writerObj, input);
    end

    % Close the video writer object
    close(writerObj);

    % 'Video.avi' will be created in the folder that contains the code.
end