clear
clf

%define some parameters
nx=100;
ny=200;
dt=1;
dx=2;                           %will have dx=dy for simplicity
Lx=nx*dx;
Ly=ny*dx;
v1=1;
v2=sqrt(2);


for i=1:nx+1
    for j=1:ny+1
        Aold(i,j)=0;
        A(i,j)=0;
    end
end

A(nx-5:nx-3,(ny/2)-1:(ny/2)+1)=200;

xx=0:dx:Lx;
yy=0:dx:Ly;
[X,Y]=meshgrid(yy,xx);
surf(X,Y,A);

%define phi
c1=(v1^2)*(dt^2)/(dx^2);
c2=(v2^2)*(dt^2)/(dx^2);

%set(gca,'nextplot','replacechildren');
%aviobj = avifile('refraction1','quality',100,'compression','Cinepak')

%h1 = figure(2)
%q = 0;

Anew=A;
for iter=1:180
    for i=2:nx
        for j=2:ny
            if i<(nx/2)+10
                Anew(i,j) = 2*A(i,j) - Aold(i,j) + c2*(A(i+1,j)+A(i-1,j)+A(i,j+1)+A(i,j-1)-4*A(i,j));
            end
            if i>=(nx/2)+10
                Anew(i,j) = 2*A(i,j) - Aold(i,j) + c1*(A(i+1,j)+A(i-1,j)+A(i,j+1)+A(i,j-1)-4*A(i,j));
            end
        end
    end
    Aold=A;
    A=Anew;
    
    if mod(iter,30)==1
        A(nx-5:nx-3,(ny/2)-1:(ny/2)+1)=200;
    end
   
    surf(X,Y,A);
    caxis([-100,100]);
    shading interp
    view(0,90)
    zlim([-200,200]);
    
    %F(j) = getframe;
    %aviobj = addframe(aviobj,F(j));
    
    %q=q+1;
    %set (gcf,'Renderer','ZBuffer');
    %if mod(iter,2)==1
    %    saveas(h1,strcat('project_movie',num2str(999+q),'.png'));
    %    fprintf('Frame %4d of %4d\n',q,180);
    %end
    rotate3d on

    pause(0.1);
end

    %aviobj = close(aviobj);