clear all
close all

%Numerical Solver for Potential Flow

N=30;
L=10;
ds = 1;
dt = .5;
X = linspace(0,L,N);
Y = linspace(0,L,N);
[x,y] = meshgrid(X,Y);

phi = zeros(N);
density = zeros(N);

density(1:4,:) = 1;

for i = 1:N
    phi(i,:) = .1*i;
end

density_new = density;
phi_new = phi;

it = 700;
for k = 1:it

    for i = 2:N-1
        for j = 2:N-1

            next_i = mod(i-1,N)+2;
            prev_i = mod(i-1,N);
            next_j = mod(j-1,N)+2;
            prev_j = mod(j-1,N);

            density_new(i,j) = density(i,j) - dt*((density(i,j)/ds^2)*(phi(next_i,j)+phi(prev_i,j)+phi(i,next_j)+phi(i,prev_j)-4*phi(i,j)) + (1/(2*ds))*((phi(next_i,j)-phi(prev_i,j))*(density(next_i,j) - density(prev_i,j)) + (phi(i,next_j) - phi(i,prev_j))*(density(i,next_j) - density(i,prev_j))));


        end
    end

    %%%%B.C.%%%%
    density_new(1,:) = density_new(2,:);
    density_new(end,:) = density_new(end-1,:);
    density_new(:,1) = density_new(:,2);
    density_new(:,end) = density_new(:,end-1);

%for n = 1:100
%    for i = 2:N-1
%        for j = 2:N-1

%            next_i = mod(i-1,N)+2;
%            prev_i = mod(i-1,N);
%            next_j = mod(j-1,N)+2;
%            prev_j = mod(j-1,N);

%            phi_new(i,j) = ((ds^2)/4)*((density_new(i,j) - density(i,j))/dt + (density(i,j)/ds^2)*(phi(next_i,j)+phi(prev_i,j)+phi(i,next_j)+phi(i,prev_j)) + (1/(2*ds))*((phi(next_i,j)-phi(prev_i,j))*(density(next_i,j) - density(prev_i,j)) + (phi(i,next_j) - phi(i,prev_j))*(density(i,next_j) - density(i,prev_j))));

%        end
%    end
%end


    density = density_new;
    phi = phi_new;

    surf(x,y,density)
    xlabel('x')
    ylabel('y')
    title(sprintf('density (iter %d)',k))
    axis([0 L 0 L -5 5])
%    view([+34.5 14]);

    if not(mod(k,4))
        saveas(gcf, ['wave_',sprintf('%04d',k),'.png']);
    end

    pause(0.001);    

end
