%% FDTD Acoustic Simulator
% Units are [m], [s], [kg], [N]
box_size=[0.1 0.1 0.1]; %simulation boundary box
gridpts=[50 50 50]; %grid points per box dim, uniformly spaced
walls_ref=0; %1 for reflecting walls (v=0), 0 for absorbing walls (P=0)

rho=1; %air density
c=314; %speed of sound

x=linspace(0, box_size(1), gridpts(1)); %convert from matrix index to physical location
dx=x(2)-x(1);
y=linspace(0, box_size(2), gridpts(2));
dy=y(2)-y(1);
z=linspace(0, box_size(3), gridpts(3));
dz=z(2)-z(1);

dt=0.5*(dx/gridpts(1))/c; %timestep

P=zeros(gridpts); %pressure matrix at (x, y, z)
Vx=zeros(gridpts); %velocity x at (x+1/2, y, z)
Vy=zeros(gridpts); %velocity y at (x, y+1/2, z)
Vz=zeros(gridpts); %velocity z at (x, y, z+1/2)

%P(gridpts(1)/2, gridpts(2)/2, gridpts(3)/2)=0.2;

%% Begin Iteration

for(i=1:20000)
    P(2, gridpts(2)/2, gridpts(3)/2)=0.2*sin(i/30); %sound wave sort of
    Vy(1:24, 20, 20:30)=0; %square tube around sound source
    Vy(1:24, 30, 20:30)=0;
    Vz(1:24, 20:30, 20)=0;
    Vz(1:24, 20:30, 30)=0;
    
    if(walls_ref==0)
        P(1,:,:)=0; %set outside P=0
    end
    Psh=circshift(P,[-1 0 0]);
    Vx=Vx-(dt/(rho*dx))*(Psh-P);

    if(walls_ref==0)
        P(:,1,:)=0; %set outside P=0
    end
    Psh=circshift(P,[0 -1 0]);
    Vy=Vy-(dt/(rho*dy))*(Psh-P);

    if(walls_ref==0)
        P(:,:,1)=0; %set outside P=0
    end
    Psh=circshift(P,[0 0 -1]);
    Vz=Vz-(dt/(rho*dz))*(Psh-P);

    if(walls_ref==1)
        Vx(end,:,:)=0; %set outside v=0
    end
    Vxs=circshift(Vx,[1 0 0]);

    if(walls_ref==1)
        Vy(:,end,:)=0; %set outside v=0
    end
    Vys=circshift(Vy,[0 1 0]);

    if(walls_ref==1)
        Vz(:,:,end)=0; %set outside v=0
    end
    Vzs=circshift(Vz,[0 0 1]);

    P=P-(rho*c^2*dt)*((Vx-Vxs)/dx+(Vy-Vys)/dy+(Vz-Vzs)/dz);
    
    if(mod(i,100)==1)
        surf(P(:,:,gridpts(3)/2));
%         hold on
%         surf(P(:,:,gridpts(3)-20)+0.5);
%         hold off
        axis([0 gridpts(1) 0 gridpts(2) -0.02 0.02]);
        title(sprintf('t=%f s',i*dt));
        pause(0.1);
        
%         frame=getframe(1);
%         img=frame2im(frame); %print(fig, '-r150', '-RGBImage');
%         [A,map]=rgb2ind(img,4);
%         if(i==1)
%             imwrite(A,map,'soundsim_Vx_center.gif','gif','LoopCount',Inf,'DelayTime',0.5);
%         else
%             imwrite(A,map,'soundsim_Vx_center.gif','gif','WriteMode','append','DelayTime',0.05);
%         end
    end
end