Tiny_FDTD_v1.m

%% Tiny_FDTD_v1.m - 7th Aug 2008 - by Nick Clark

close all;

%% Some user modifiable parameters
xDim = 80;      %set number of nodes in x direction
yDim = 80;      %set number of nodes in y direction
c = 343;        %speed of sound m/s
rho0 = 1.21;    %density of air kg/m^3
dx = 0.1;       %distance between nodes

t0 = 100;       %index time of gaussian pulse peak
width = 10;     %peakyness of gaussian pulse

%% Initialise the p and u fields and set coefficients + a Courrant compliant dt
p = zeros(yDim,xDim);       %pressure scalar field initialisation
ux = zeros(yDim,xDim-1);    %velocity vector field initialisation
uy = zeros(yDim-1,xDim);    %velocity vector field initialisation

dt = dx/(sqrt(ndims(p))*c); %set time step small enough to satisfy courrant condition
coP = (dt/dx) * rho0 * c^2; %pressure scalar difference co-efficient
coU = (dt/dx) / rho0;       %velocity vector difference co-efficient

%% FDTD loop - comment out plot lines for speed
for n = 1:400;
    uy = uy + coU.*( p(1:end-1, :) - p(2:end, :) );              %calculate velocity in y plane
    ux = ux + coU.*( p(:, 1:end-1) - p(:, 2:end) );              %calculate velocity in x plane
    p = p + coP.*( [zeros(1,xDim); uy] - [uy; zeros(1,xDim)] ...
                 + [zeros(yDim,1) ux] -  [ux zeros(yDim,1)] );   %calculate pressure from xy velocity
    p(ceil(yDim/2),ceil(xDim/2)) = exp(-.5 * ((n-t0)/width).^2); %insert hard source
          
    surf(p); shading interp; lighting phong; colormap hot; axis off; zlim([0 1]);
    set(gcf,'Color', [0 0 0], 'Number', 'off', 'Name', sprintf('Tiny FDTD, step = %i', n));
    title(sprintf('FDTD Time = %.2f msec',n*dt*1e3),'Color',[1 0 0],'FontSize', 22); drawnow;
end;