% Number of cells and timesteps
nx = 50;
ny = 25;
nt = 5;

% Size of total domain
width = 100;
height = 100;

% Size of each cell
dx = width / nx;
dy = height / ny;

% Wave speed
c = 1.0;

% Initial heat distribution
u0 = zeros(ny, nx);
u0(12, 25) = 1;
u0(11:13, 24) = 0.5;
u0(11:13, 26) = 0.5;
u0(11, 25) = 0.5;
u0(13, 25) = 0.5;
u1 = u0;
u2 = u0;

% Generate x and y coordinates for each cell
x = linspace(0.5*dx, width-0.5*dx, nx);
y = linspace(0.5*dy, height-0.5*dy, ny);
[x, y] = meshgrid(x, y);

% Maximum size of timestep (according to CFL)
cfl = 0.8;
dt = cfl*min(dx*dx/(2*c), dy*dy/(2*c));

% Plotting help
figure('units','normalized','outerposition',[0 0 1 1]);
simulated_data = surf(x, y, u2);
zlim([-max(max(u2)), max(max(u2))]);
legend('Pressure');
title('Linear wave equation in 2D');

% Loop over time for nt timesteps
for k=0:nt
    
    % Loop over all the internal cells
    for j=2:ny-1
        for i=2:nx-1
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Dette er linjen du må endre på %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            u2(j, i) = u1(j, i);
        end
    end
    
    % Set reflective boundary conditions
    u2(1, 1:nx) = u2(2, 1:nx);
    u2(end, 1:nx) = u2(end-1, 1:nx);
    u2(1:ny, 1) = u2(1:ny, 2);
    u2(1:ny, end) = u2(1:ny, end-1);
    
    % Rotate / swap the data
    u0 = u1;
    u1 = u2;
    
    % Update plot
    simulated_data.ZData = u2;
    drawnow;
    pause(0.01);
end



